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SECTION  I 


INTRODUCTION 

1,1  Objectives  and  Rationale 

This  volume  contains  the  analytical  development  of  the  equations 
required  for  mechanization  of  a strapdown  inertial  navigation  system 
simulation  program. 

One  objective  of  this  volume  is  to  help  clarify  the  meanings  of 
the  relationships  between  the  multiply-defined  coordinate  frames  which 
occur  throughout  the  program.  These  multiply-defined  frames  do  not 
materially  complicate  the  coding  but  can  impede  analytical  insight  un- 
less fully  understood.  This  effort  is  concentrated  in  the  Appendices 
(mainly  Appendix  B)  and  the  results  are  simply  introduced  in  the  main 
body. 


Another  objective  is  to  provide  the  necessary  analytical  foundation 
for  the  equations  which  often  merely  appeared  in  earlier  documentation. 
Program  structure  was  unchanged  but  several  modules  were  reworked t the 
laser  gyro  and  compensation  are  completely  new?  the  navigation  and  velo- 
city-attitude algorithm  are  upgraded  per  (5). 

The  approach  to  providing  the  analytical  development  of  the  equations 
actually  coded  started  with  inverse  Fortran  transformation  of  the  code  - 
to  obtain  the  difference  equations  mechanized  - followed  by  rederivation 
of  these  equations  and  corrections  where  necessary. 

1.2  Suxttnary 

Starting  with  discussion  of  the  need  filled  by  the  INSS  and  the 
general  structure  of  the  program,  the  remainder  of  the  volume  derives, 
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n 

develops,  and  the  presents  revelant  algorithms,  generally  on  a module- 

A 

by-module  basis. 


Section  3 describes  the  interim  trajectory  module-a  self-contained 
driver  which  does  not  permit  aircraft  maneuvers.  Although,  the  pre- 
decessor of  this  module  was  used  almost  exclusively  in  the  earlier  version 
of  program,  it  was  not  documented  in  the  reports.  This  time  the  trajectory 
module,  which  interfaces  with  PROFGEN,  is  omitted  since  it  and  a modified 
PROPGEN  program  are  currently  available  at  AFAL,. 


Section  4 develops  and  presents  the  equations  and  algorithms  nec- 
essary to  simulate  directly  the  random  linear  and  angular  motion  of  the 
vehicle  defined  by  six  displacement  power  spectral  densities.  The 
resultant  random  motion  is  superimposed  on  the  trajectory  output  and 
passed  on  to  the  inertial  component  modules  but  not  to  the  evaluation 
module,  hence  vibration  appears  as  navigation  errors. 

Section  5 develops  and  presents  the  models  for  a triad  of  single* 

degree-of -freedom,  floated,  rate  integrating  qyroa  (operated  in  a rate 

mode)  or  ring  laser  gyros  and  the  related  compensation,  which  latter  is 

performed  in  the  simulated  navigation  computer.  The  dynamics  of  the  SOOF 

gyros  are  derived  from  first  principles,  then  reduced  to  algorithmic  form, 

where  the  user  is  given  the  choice  of  three  dynamic  models  varying  in 

effective  bandwidth.  The  treatment  of  component  alignment,  is  rather 

complex,  and  despite  the  explanation  (in  Appendix  P) , must  be  approached 

with  caution.  The  laser  gyros  and  compensation  models  are  completely  new 
&**■ 

and  -th*  baaed  mainly  on  test  data.  They  are  vastly  simplified  relative  to 
their  predecessors  in  the  earlier  versions  of  the  program.  The  user  tnay 
select  either  type  of  gyro  by  appropriate  module  substitution. 

Section  6 develops  and  presents  the  models  for  a triad  of  strapdown, 
pendulous,  single  axis,  viscous  damped  accelerometers  with  integrated, 
digital  outputs.  The  linear  acceleration  resulting  from  .locating  the 
accelerometers  some  distance  from  the  center  of  rotation  of  the  vehicle 
is  included,  in  addition  to  the  errors  sources  affecting  SDOF,  floated 
accelerometers.  Three  dynamic  models  may  be  selected.  The  compensation 
scheme,  which  is  included,  is  more  complete  than  that  appearing  in  most 
"operational**  mechanisation. 
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Section  7 presents  the  velocity-attitude  computations  for  a strap- 
down  inertial  navigator.  Since  the  relevant  equations  and  algorithms 
were  developed  in  (3)  and  (5) , there  was  no  necessity  to  develop  them 
from  first  principles.  Updating  of  the  alignment  qviaternion  or  direction 
cosine  matrix  (d.c.m)  by  a first,  second,  or  third  order  algorithm  is 
permitted,  with  ( ortho)  normalization  at  a user-specified  frequency. 
Regardless  of  the  form  of  update  employed,  a d.c.m.  is  formed,  since  it 
is  required  in  attitude  extraction.  A separate,  mid-computation  cycle 
d.c.m.  is  computed  for  incremental  velocity  transformation. 

Section  d presents  the  navigation  computations  for  a strapdown  INS, 
wherein  incremental  velocity  is  summed  in  an  inertial  frame,  then  trans- 
formed to  a LVWA  frame  for  integration.  The  navigation  algorithm  cori'es- 
ponds  closely  to  the  upgraded  algorithm  of  (5),  except  for  the  incremental 
velocity  summation  in  the  inertial  rather  than  the  LVWA  frame.  Attitude 
is  extracted  from  the  appropriate  matrix  product.  Third  order  vertical 
channel  damping  is  employed  using  the  simulated  barometric  altimeter  input 
as  a reference. 

Section  9 briefly  describes  the  error  evaluation  routine. 

Appendix  A,  extracted  from  (S)  defines  single  axis,  rotations  and 
the  general  rotation  matrix. 

Appendix  B discusses  and  describes  the  simplest  set  of  coordinate 
frames  and  transformations  which  could  be  employed  in  tin*  simulation  of 
a strapdown  INS  with  a LVWA  computational  frame,  and  transformations 
actually  employed. 

Appendix  C presents  the  rotation  matrix  in  toms  of  the  direction 
cosines  of  the  axis  of  rotation  and  the  Angle  of  rotation  and  shows  the 
relationship  between  a d.c.m  and  a unit  quaternion. 

Appendix  0 discusses  the  inertial  component  alignment  matrices 
and  indicates  how  the  required  117  entries  describing  component  alignment 
may  be  generated  for  the  component  models  and  compensation. 
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SECTION  2 
OVERVIEW  OP  INSS 


2.1  The  Evolution  of  the  Strapdown  Inertial  Navigation  System 

2.1.1  Historical  Prologue 

Twenty  years  ago,  inertial  navigation  and  guidance  systems  were  in 
their  infancy. 

The  inertial  measuring  unit  (IMU)  contained  three  single-degree- 
of- freedom  rate  integrating  gyros  and  two  or  three  floated  pendulous 
accelerometers  or  pendulous  integrating  gyro-accelerometers  (PIGA) . 
Between  the  inertial  sensors  and  the  IMU  case  were  three  to  five  gim- 
bals, each  with  its  own  gimbal  torque  motor  (usually  geared,  two-phase) 
and  one  or  more  resolvers  and/or  synchros. 

The  navigation  con\puter  was  inevitably  analog,  generally  with 
electromechanical  integrating  servos,  employing  motors  with  tachometer- 
generator  feedback.  The  platform  resolvers  were  frequently  part  of  the 
analog  computer.  Vacuum  tubes  were  being  replaced  by  transistors;  the 
silicon  power  transistor  was  brand  new.  Coordinate  transformations  (if 
required)  employed  resolver  chains.  The  first  airborne,  transistorized, 
digital  computer  - with  a digital  differential  analyser -was  in  the  offing. 

The  mechanization  of  choice,  for  terrestrial  (aircraft)  applica- 
tions, was  the  "analytic"  or  local-level,  north-pointing  system,  where- 
in three  of  the  gimbal  synchros  provided  direct  indication  of  the  air- 
craft attitude  - roll , pitch  and  heading  (in  the  five  gimbal  "geometric" 
system,  synchros  or  encoders  also  provided  direct  outputs  of  latitude 
and  longitude). 

Two  other  possible  mechanizations  were  of  largely  academic  inter- 
est at  this  time,  the  "space  stable"  (or  inertially-stabilizod)  system 
and  the  "gimballoss"  (or  strapdown)  system.  Both  of  those  mechaniza- 
tions were  awaiting  tho  advent  of  smaller,  faster,  more  powerful,  air- 
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borne  digital  computers  and  the  development  of  pulse-rebalanced  gyros 
and  accelerometers. 

2. 2 The  Dawn  of  the  Strapdown  System 

With  the  advent  of  the  sixties,  active  strapdown  system  development 
began  in  earnest.  Some  of  the  advantages  of  the  gimballed  mechaniza- 
tions soon  came  to  be  more  widely  appreciated,  namely  the  isolation  of 
the  inertial  components  from  the  angular  motion  environment,  by  the 
gimbal  servos  at  low  frequencies  and  by  the  inertia  of  the  stable 
element  at  higher  frequencies.  Generally,  further  isolation  of  the 
inertial  components  from  angular  and  linear  vibration  environment  was 
provided  by  platform  isolators , with  natural  frequencies  on  the  order 
of  20  to  50  Hz. 

Most  of  the  structural  modes  of  the  aircraft  are  below  20  Hz,  so 
a strapdown  system,  even  with  isolators,  would  be  subjected  to  sizable, 
coupled,  angular  and  linear  vibration  inputs  at  low  frequencies,  as 
well  as  the  gross  vehicle  motion  due  to  maneuvers.  Some  of  the  result- 
ing errors  are  due  to  the  kinematics  of  the  motion  experienced  by  the 
sensors,  as  in  "coning"  and  "anisoinertia" , and  will  affect  even  per- 
fect inertial  sensors,  while  others  are  the  results  of  the  interaction 
of  the  environment  and  the  dynamics  of  the  particular  sensor  and  its 
rebalance  loop.  Only  in  very  special  cases  are  the  effects  of  the 
above  types  of  forcing  functions  on  strapdown  inertial  navigation  sys- 
tem performance  analytically  tractable.  Thus  the  available  methods  of 
assessing  strapdown  inertial  component  or  system  performance,  short  of 
actual  flight,  are  reduced  to  environmental  test  and  simulation. 

Strapdown  systems  also  impose  rather  severe  computational  require- 
ments on  their  navigation  computers.  These  arise  from  the  necessity  of 
transforming  the  incremental  velocity  outputs  of  the  accelerometer  from 
the  vehicle  or  body  frame  to  an  inertial  or  earth-referenced  computa- 
tional frame  in  which  the  position  and  velocity  of  the  vehicle  are 
computed. 
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Since  this  transformation  includes  the  total  vehicle  angular 
motion  due  to  both  maneuvers  and  vibration  (up  to  several  hundreds  of 
degrees  per  second) , it  has  to  be  computed  very  frequently  and  very 
accurately.  Computational  errors  in  strapdown  systems  could  easily 
equal  or  exceed  the  entire  system  error  budget.  Here  again,  the  analy- 
tical assessment  of  computational  error  propagation  in  strapdown  algo- 
rithms generally  included  only  those  special  cases  which  were  analy- 
tically tractable. 

These  analytical  results  led  to  the  re-emergence  of  the  quaternion 
(from  the  mists  of  antiquity)  as  the  preferred  means  of  generating  the 
transformation  from  the  body  (inertial  sensor)  frame  to  the  computation- 
al frame.  Very  recent  indications  are  that  at  least  one  strapdown  sys- 
tem manufacturer  is  returning  to  the  "old"  direction  cosine  or  rotation 
matrix  directly. 

2 . 3 The  Present  Situation 

The  first  strapdown  system  was  flown  about  ten  years  ago.  Today, 
there  are  a number  of  strapdown  systems  in  varying  stages  of  develop- 
ment and/or  test.  Flight  test  results  range  from  low  accuracy  to  mod- 
erately high  accuracy.  With  regard  to  the  inertial  conponents,  all  the 
accelerometers  used  are  of  the  pendulous,  single-axis  variety,  gener- 
ally with  analog  loops  and  external  analog-to-digital  conversion.  A 
full  range  of  gyro  types  are  found  ranging  from  two-degree-of-freedom, 
free  gyros  (ESG's)  and  two-degree-of-freedom,  dry-tuned  gyros  to  single- 
degree-of-freedom  floated,  rate-integrating  gyros  and  laser  gyros,  with 
a variety  of  gyro  signal  digitization  methods.  Only  the  laser  gyro  has 
an  inherently  digital  output. 

The  digital  computers  for  the  strapdown  IMU's  are  changing  so  rapid- 
ly that  no  general  characterization  is  likely  to  be  correct  for  long, 
beyond  the  fact  that  the  word  length  is  ususally  16  bits  and  the  cycle 
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time  is  less  than  one  microsecond.  They  are  small,  fast,  powerful,  and 
cheap,  compared  with  their  raini-based  predecessors  of  a few  years  ago. 

2 . 4 Justification  for  the  Inertial  Navigation  Strapdown  Simulator  (INSS) 

With  the  proliferation  of  existing  and  projected  strapdown  INSS's, 
the  need  for  a flexible,  powerful  tool  for  the  evaluation  of  the  pro- 
jected performance  of  candidate  strapdown  systems  in  their  operational 
environment  became  evident. 

It  was  in  an  attempt  to  fulfill  this  need  that  the  INSS  was 
generated. 


2.5  Characteristics  of  the  INSS,  or  What  Does  the  Simulator  Do? 

The  simulator  performs  the  following  functions: 

a)  Provides  a flight  profile  for  a point-mass  aircraft  executing 
coordinated  maneuvers.  Profile  includes  position,  velocity, 
acceleration  and  attitude  of  the  aircraft. 

b)  Superimposes  on  the  flight  profile,  the  angular  and  linear 
motion  of  the  aircraft  in  response  to  random  aerodynamic 
forces  such  ~s  mist*  and  turbulence. 

o)  Operates  on  the  specific  force  and  angular  velocity  to  produce 
the  ideal  body  frame  values  of  these  quantities  and  the  corres- 
ponding angular  acceleration. 

d*  Integrates  the  gyro  and  accelerometer  equations  of  motion,  of 
the  form  specJ*ied,  after  accounting  for  the  displacements  of 
the  inert! --I  components  from  the  "center"  of  the  point-mass 
vehicle,  and  their  orientations  with  respect  to  the  body  frame. 
Incorporates  effects  of  component  parameters  and  environment 
on  component  outputs. 
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e)  Reads  true  altitude  and  perturbs  the  same  according  to  the  al- 
timeter model. 

f)  Transmits  sensor  data  from  simulated  IMU  to  simulated  naviga- 
tion computer,  and  resets  sums. 

g)  Performs  accelerome ter  compensation  function  (this  is  essen- 
tially the  inverse  of  the  accelerometer  model,  with  errors  and 
omissions) . 

h)  Using  compensated  accelerometer  outputs  (Av's)  and  gyro  out- 
puts (A0's)  it  performs  the  gyro  compensation  (again,  this  is 

essentially  the  inverse  of  the  gyro  model,  with  errors  arid 
omissions) . 

i)  Updates  the  body- to- inertial-frame  transformation  (direction 
cosine  matrix,  DCM)  and  transforms  the  incremental  velocities 
to  the  inertial  frame  (either  a DCM  or  a quaternion  update 
may  be  employed) . 

j)  Transforms  the  incremental  velocities  to  the  local  vertical 
wander  azimuth  computational  frame  and  aomputes  local  vertical 
position  and  velocity  (incorporating  the  barometric  altimeter 
for  vertical  danping)  and  attitude. 

k)  As  required,  the  navigator  outputs  are  differenced  with  the 
flight  profile  values  and  the  resulting  errors  are  printed 
out  and/or  plotted. 

These  functions  are  shown  pictorially  in  Figure  2-1, 
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INERTIAL  NAVIGATION  STRAPOQWN  SIMULATION  PROGRAM 


Figure  2-1.  INSS  FUNCTIONAL  BLOCK  DIAGRAM 


SECTION  3 


TRAJECTORY  GENERATION 


The  INSS  may  employ  any  trajectory  module  which  is  compatible  with 
the  conventions  used  in  the  remainder  of  the  modules. 

The  trajectory  module  which  will  normally  be  used  to  investigate 
the  effects  of  vehicle  dynamics  encountered  in  tactical  aircraft  is 
available  from  the  Air  Force  Avionics  Laboratory  (RWA-3).  This  module 
accepts  the  outputs  of  an  independent,  AFAL-deve loped , flight  profile 
generator  program,  PROFGEN (2) , which  has  recently  been  modified  to 
simulate  more  closely  the  dynamics  of  a point-mass  aircraft. 

PROFGEN  outputs  will  include  the  geodetic  position  (and  wander 
angles  if  appropriate)  earth-relative  velocity,  and  attitude  (roll, 
pitch,  and  yaw  or  heading) , in  addition  to  specific  force  (or  incre- 
mental velocity)  in  any  of  several  preselected  coordinate  frames.  The 
earth  and  gravity  models  employed  in  PROFGEN  are  based  on  the  WGS72 
Ellipsoidal  Earth  Model  (1)  and  hence  are  compatible  with  the  INSS 
modules. 

The  function  of  the  INSS  trajectory  module,  is  this  case,  is  largely 
a matter  of  accepting  the  PROFGEN  outputs  and  performing  such  operations 
(as  initialization,  summation,  transformation,  etc)  as  are  necessary 
before  passing  the  data  on  to  the  remaining  INSS  modules.  The  above 
trajectory  modulo  and  the  flight  profile  generator  program,  PROFGEN,  are 
not  described  further  in  this  report.  Instead,  the  interim  trajectory 
module  (which  provides  the  same  output  data)  is  described  in  the  remainder 
of  this  section. 
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The  trajectory  generation  module  (TRJ)  allows  the  user  to  simulate 
elementary  flight  profiles  by  specifying  a simple  set  of  flight  control 
parameters.  The  module  TRJ  was  designed  as  a program  driver  for  software 
development  and  debugging;  therefore,  it  has  very  limited  capabilities. 

It  is  restricted  to  constant  velocity  profiles  at  fixed  aircraft  headings. 
In  addition,  there  is  no  provision  for  roll  or  pitch  motion. 

The  user  specifies  the  initial  position,  wander  angle,  and  attitude 
of  the  aircraft,  plus  the  east,  north,  up  components  of  velocity  relative 
to  the  earth. 


V 


These  velocity  components  are  constant  throughout  the  simulated  flight. 
Given  the  initial  conditions,  the  program  computes  the  current  values  of: 


• Position 

• Velocity  (constant) 

• Heading  (constant) 

• Roll,  pitch,  yaw  (R  <=  const,  P » const) 

• Azimuth  wander  angle 

• Integrated  specific  force  in  body  coordinates 

• Specific  force  in  body  coordinates 

• Body  angular  velocity  in  body  coordinates 


L,  l,  h 


R,  P,  V 


a 


<ob 

-ib 


The  specific  force  and  angular  velocity  are  sent  to  the  random 
environment  module  (ENV),  She  altitude  and  vertical  velocity  are  used 
in  the  altimeter  module  (ALT!) . She  evaluation  module  (EVL)  compares 
the  values  of  position,  velocity  and  attitude  from  the  trajectory  module 
with  the  values  computed  by  the  navigation  algorithms. 


The  trajectory  nodule  reads  two  files  of  input  data  when  it  is 
initialized.  During  the  initialization  pass  it  also  generates  an  output 
file  containing  geodetic  constants  and  initialization  data  which  is 
available  to  the  other  program  modules.  The  initialization  of  the 
trajectory  module  is  summarized  in  Section  3*1.  The  normal  module 
computations  are  described  in  Section  3.2. 

3.1  Trajectory  Module  Initialization 

The  following  operations  are  performed  on  the  initialization  pass 
through  the  trajectory  module. 

1.  Read  the  two  initialization  data  files,  XFXLE  and  PFILE,  and 
convert  the  data  into  internal  program  units. 


XFXLE  INPUT  DATA 

Variable 

Description 

Units 

DT 

module  time  step 

s 

PRNTSW 

print  control  switch 

- 

QUTSW 

not  used 

- 

XFXLE 

unit  number  of  print  file 

«■* 

ILAT 

geodetic  latitude 

deg 

I LON 

geodetic  longitude 

deg 

I ALT 

altitude 

ft 

XVEL  (1) 

ground  velocity-east 

ft/s 

XVEL  (2) 

ground  velocity-north 

ft/s 
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Variable 

Description 

Units 

IVEL  (3) 

ground  velocity-  up 

ft/s 

IPITCH 

pitch  angle 

deg 

IYAW 

yaw  angle 

deg 

IROLL 

roll  angle 

deg 

MODPDT 

module  print  interval 

s 

IWANDER 

initial  azimuth  wander  angle 

deg 

PFILE  INPUT  DATA 


Variable 

Description 

Units 

WE 

earth  rotation  rate 

rad/s 

RE 

equatorial  earth  radius 

ft 

G 

nominal  gravity 

ft/s2 

PRNTDT 

general  print  interval 

s 

2.  Compute  the  specific  force  in  body  coordinates.  The  spocific 
force  in  geographic,  i.o.,  (e.  n,  u) , coordinates  is 


^ . n . 

Z a * <%n  + 2ujie)  x vj  - a 


t 


where  V is  the  velocity  relative  to  the  earth  and  £ is  the  gravity 

vector,  Tho  trajectory  module  is  restricted  to  constant  velocity,  so 
. i 

V is  zero.  Tho  program  also  assumes  the  velocity  V is  zero  during  the 

“**  G *"*6 

initialization  pass  and  computes 


- 


j>  a 
a c a 
n 
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This  is  correct  only  if  the  earth  relative  velocity  is  zero. 
Fortunately,  the  accelerometer  and  gyro  modules  do  not  use  the  initial 
value  of  specific  force,  so  the  approximation  is  immaterial. 

3.  Compute  the  angular  velocity  of  the  body  with  respect  to 
inertial  space  in  body  coordinates.  Since  the  orientation  of  the  body 
is  fixed  relative  to  geographic  coordinates,  the  angular  velocity  of 
the  body  wrt  inertial  space  is  equal  to  the  angular  velocity  of  the 
geographic  frame  (the  2-frame)  wrt  inertial  space. 

Hib  m HU 
" Hie+ HU 


The  angular  velocity  expressed  in  inertial  coordinates  is 


"o  ’ 

'-(V/r)  cos  l ' 
n m 

wie 

V /(r  cos  L) 
e P 

0 

■ « 

(V  /r  ) sin  £ 

U n m j 

where  r and  r are  the  meridional  and  prime  vertical  radii  of  curvature 
m p 

The  angular  velocity  is  then  transformed  to  body  coordinates. 


J>  i 

ciAb 


4.  Compute  the  aircraft  heading  from  the  specified  azimuth  wander 
and  yaw  angles.  The  heading  remains  constant  throughout  the  flight 
trajectory. 


H ■>  ¥ - a 


where  ¥ is  the  specified  initial  yaw  angle  and  a is  the  initial  wander 
angle. 

5.  Compute  the  time  rate  of  change  of  the  roll,  pitch  and  yaw 
angles.  Singe  the  body  orientation  is  fixed  relative  to  geographic 
coordinates, 
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R 


0 


(roll  rate) 


* 

P = 0 (pitch  rate) 

Ve 

Y = a tan  L (yaw  rate)  . 

P 

6.  Generate  and  write  the  output  data  file  PF1LE.  The  output 
data  consists  of  initialisation  constants  from  the  IFILE  and  PFILE  input 
data  files,  plus  the  initial  value  of  specific  force  which  was  computed 
during  the  initialization  pass.  The  output  PFILE  is  available  to  the 
other  program  modules.  Although  the  roll,  pitch  and  yaw  rates  are 
placed  in  the  PFILE,  they  are  not  used  by  any  of  modules.  All  variables 
in  the  PFILE  are  in  internal  program  units. 

PFILE  OUTPUT  DATA 


Variable 

Description 

Units 

WE 

earth  rotation  rate 

rad/s 

RE 

equatorial  earth  radius 

ft 

1 

nominal  gravity 

ft/s2 

PRNl’DT 

general  print  interval 

s 

LAT 

geodetic  latitude 

rad 

LON 

geodetic  longitude 

rad 

WANDER 

initial  azimuth  wander  angle  rad 

ALT 

altitude 

ft 

ROLL 

roll  anglo 

rad 

PITCH 

pitch  angle 

rad 

YAW 

yaw  angle 

rad 

ROOT 

roll  rate 

rad/s 

poor 

pitch  rate 

rad/s 

YOGT 

yaw  rate 

rad/s 

VEL(l) 

ground  velocity-east 

ft/s 

VEL(2) 

ground  velocity-north 

ft/s 

VEL{3) 

ground  velocity-up 

ft/s 

ABU) 

1 K 1 

ft/s2 

ABU) 

t specific  force  , 

ft/s2 

AB(3) 

( ' I 

ft/s2 
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3.1.1  Position  Equations 


The  program  is  restricted  to  constant  velocity  with  respect  to 
the  earth  in  geographic  coordinates.  The  module  mechanizes  a conven- 
tional set  of  equations  to  compute  latitude,  longitude  and  altitude 

from  the  specified  velocity  (V  • V , V ).  The  differential  equations 

6 n u 

are: 


w 


ee 


-V  /r 
n'  m 


U) 


en 


VrP 


to  * (V  /r  ) tan  L 
eu  e p 


i * to/  cos  L 
en 

* • vu  ‘3-11 

where  rw  and  r^  are  tho  meridional  (northerly)  and  prime  vertical 
(easterly)  radii  of  curvature  respectively. 


reU  - e ) 


(1  - c2  sin2  U 3/2 


tl  - e2  sin2  UW 


geodetic  latitude 
geodetic  longitude 
altitude 


h 


+ h 


* {vv\} 

■ 

» equatorial  earth  radius  (20325640  ft) 
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e 


eccentricity 


e2  = e(2  - e)  » 0.006694317778 
e = 1/298.26  = ellipticity 

The  differential  equations  are  solved  using  rectangular  integrations 


n + 1 


n + 1 


n + 1 


= L - Ul  At 
n ee 


= 2 + (m  /cos  L)At 

n en 


a h + V At 
n u 
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3.1.2  Attitude 

« • 

There  is  no  roll  or  pitch  motion;  i.e.,  R and  P equal  zero.  The 
heading  is  also  constant.  The  yaw  angle  in  the  program  is  defined  as 


V • a - ii 


so,  the  yaw. rate  is  equal  to  the  rate  of  change  of  the  azimuth  wander 
angle. 


-w 


cu 


Rectangular  integration  is  used  to  compute  the  yaw  and  azimuth 
wander  angles. 


a 

n 


+ 1 


a - w At 
n ou 


¥ 

n 


♦ 1 


V 

n 


w At 

QU 
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The  roll,  pitch,  yaw  and  azimuth  wander  angles  are  used  to  com- 

i>  b 

pute  the  body  attitude  matrices  Cj,  and  CL.  Hie  transformations  utilize 
the  local  level  azimuth  wander  frame  (c). 
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< 


(3-4) 

<4 

(3-5) 

The  C transformation  is  a function  of  the  roll,  pitch,  and  yaw  angles, 
c 

Cb  = X (II)  X (R)  Y (P)  Z (Y)  X (II)  Z (~) 
c * 


The  transformation  from  (£)  to  the  azimuth  wander  frame  is  a simple 
rotation  through  the  azimuth  wander  angle, 

4 =■  Z.(0) 

The  details  of  each  of  the  transformations  is  explained  in  Appendix  B. 

Notice  that  while  the  program  uses  the  azimuth  wander  frame  (c) , 
this  coordinate  system  is  not  needed.  It  would  have  been  simpler  to 
transform  directly  to  the  body  frame  using  the  heading  angle. 

Cb  = C(R)C(P}C(H) 


3.1.3  Angular  Velocity  of  the  Body 

The  average  angular  velocity  of  the  body  over  the  simulation  time 
step  At  is  computed  from  the  Euler  rotation  of  the  body  over  the  time 
interval. 


Cb{n)C^(n 


- 1}  « M(At) 


-F(0b ) 

W Q 


« I 


- F(j^) 


where 
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The  incremental  Euler  rotation  is  obtained  from  the  off-diagonal  elements 
of  the  matrix  M;  viz0 , 


\ 
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The  average  angular  velocity  in  body  coordinates  is 

uk  = 9b/At 
—lb  — 
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3.1.4  Specific  Force 

The  specific  force  in  geographic  coordinates  is  computed  from 
the  specified  velocity.  The  specific  force  is  then  transformed  to  the 
body  frame.  The  general  expression  for  the  specific  force  in  (e,  n,  u) 
coordinates  is 


a* 


V' 


(u) 


,enu 


2^E> 


x V 


The  program  computes  the  integral  of  specific  force  over  the  simulation 
time  step  At.  The  integrated  specific  force  is  made  available  to  the 
evaluation  module  (EVL)  for  comparison  purposes. 


dV2  = AV2  + [ (U)  . + 2U>.  ) x V2  - g2]At 
— — e — e x.  — le  — « 


The  scalar  expressions  are 


- (U)  + 2u>,  sin  L)V  + u>  + 20).  cos  L)V 

eu  ie  n en  re  u 

(0)  + 2u).  sin  L)  V - u)  V + g 

eu  ie  e ee  u n 


„ - ( 0)  + 2 0).  cos  L)  V + 

L en  ie  e 


0)  V + g 
ee  n u 


At 
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where  the  components  of  angular  velocity  with  respect  to  the  earth 

(<ee*  %'  “eu5  are  given  in  Eg‘  (3-1).  The  WGS72  coefficients  for 
the  gravity  vector  are 

9n  “ 1.63  * 10  8 h sin  L cos  L 

9U  - -(32.0877057  + 0.16939081  sin2  L + 7.52810  x lO-4  sin4  L) 

x U - (9.6227  X 10  8 - 6.4089  x lo”10  sin2  L)  h + 6.8512 
x 10"i5h2] 

The  change  in  velocity  AV^  is  zero  because  the  specified  velocity 
is  constant.  The  specific  force  is  then  computed  from  dv*  and  trans- 
formed to  body  coordinates. 


(3-9) 
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SECTION  4 


RANDOM  ENVIRONMENT 


The  environment  module  . (ENV)  simulates  angular  and  linear  vehicle 
motion  generated  by  unpredictable  sources  such  as  wind  buffeting  during 
ground  alignment  and  air  turbulence  during  flight.  The  current  version 
of  the  program  is  restrictive  because  it  allows  only  one  representation 
of  random  motion  for  all  segments  of  the  simulated  flight  profile.  This 
deficiency  should  be  corrected  because  the  environment  during  ground 
alignment  is  drastically  different,  for  example,  from  the  inflight  en- 
vironment of  high-performance  aircraft.* 

Power  spectral  density  plots  are  usually  used  to  describe  random 
aircraft  motions.  The  shape  of  the  PSDs  is  usually  quite  complex.  For 
example,  the  PSD  in  Figure  4-1  shows  the  motions  measured  in  the  B-l 
radome  - forward  avionics  bay. 

The  environment  module  accepts  PSD  characterizations  for  each 
component  {roll,  pitch,  yaw)  of  linear  and  angular  motion.  The  module 
converts  the  PSD  representations  into  linear  and  angular  random  motions 
as  functions  of  time.  In  designing  the  program  it  was  assumed  each  of 
the  six  components  of  random  motion  were  uncorrelated.  While  this 
assumption  is  certainly  not  valid,  it  is  convenient  because  cross- 
correlation information  is  normally  not  available.  The  program  accepts 
PSDs  at  the  displacement  level;  i.e. 

2 

• Angular  displacement  PSD  - (rad  /Hz) 

2 

• Linear  displacement  PSD  - (ft  /Hz) 


♦Alternatively,  the  capability  to  run  the  program  in  segments  e.g 
ground  alignment,  cruise,  etc.,  could  be  added. 


The  position  level  PSDs  must  be  constructed  if  the  aircraft  motion  is 
specified  at  the  velocity  or  acceleration  level. 

To  simulate  the  random  motion  it  is  necessary  to  specify  the 
following  items: 

a.  Descriptions  of  the  Six  PSDs 


Each  PSD  is  approximated  by  specifying  parameters  which  describe 
the  PSD  in  each  regions  of  peak-power  density.  Three  parameters 
are  needed  for  each  region  of  peak  density  (see  Figure  4-1) . 

2 

• Ap  - peak  PSD  amplitude  (units  /Hz) 

• BW  - half  power  bandwidth  (rad/s) 

• w - frequency  at  peak  (rad/s) 

The  number  of  peaks  must  also  be  specified.  The  program  allows 
a maximum  of  five  peaks  for  each  PSD. 

* 

Wind-Gust  Intensities 

The  PSD  amplitudes  specified  above  represent  vehicle  motion 
when  the  rms  gust  intensity  along  each  axis  is  1 foot/second. 

The  variances  of  the  desired  intensity  levels  must  be  speci- 
fied. 

• Lateral  (ft/s)2 

2 

• Longitudinal  (ft/s) 

• Normal  (yaw)  (ft/a)2 

simulation  Time  Ste 


The  simulation  time  step  must  be  chosen  so  the  simulation 
frequency  is  approximately  10  times  the  highest  frequency 
of  interest.  The  program  will  generate  random  motions  which 
approach  discrete  white  noise  if  the  simulation  frequency  is 
too  low.  The  highest  frequency  should  be  obtained  from  angular 
velocity  PSDs  and  linear  acceleration  PSDs  because  the  position 
level  plots  attenuate  the  high-frequency  components. 


4.1  Random-Motion  Computations 

Figure  4-2  shows  the  operations  performed  in  the  random  motion 
module.  The  linear  displacement  PSDs  are  converted  into  samples  of 
random  velocity  in  body  coordinates.  The  velocity  samples  are  then 
differenced  to  form  samples  of  angular  accelerations.  The  angular 
displacement  PSDs  are  converted  into  samples  of  angular  velocity  in 
body  coordinates.  The  average  angular  velocity  <$u  over  two  simulation 
cycles  is  used  to  compute  the  transformation  from  nonvibrating  to 
vibrating  body  coordinates. 

The  total  linear  acceleration  is  obtained  by  adding  the  random 
acceleration  6ab  to  the  acceleration  from  the  flight-profile  generator 
after  it  has  been  transformed  to  vibrating  body  coordinates.  The  total 
angular  velocity  is  obtained  by  adding  the  random  angular  velocity  6ub 
to  the  angular  velocity  from  the  flight-path  generator.  The  derivative 
of  angular  velocity  is  formed  by  differencing  successive  samples  of  total 
angular  rate. 

4.1.1  Random-Motion  Generator 

The  algorithms  for  the  linear  and  angular  random  motion  generators 
are  identical  - only  the  units  are  different.  Each  displacement  PSD  is 
converted  into  random  displacement  motion  <SX  samples. 

6vn  “ <6xn  " 6Vl)/At 

Displacements  are  simulated  rather  than  velocity  or  acceleration  to 
guarantee  that  random  displacements  will  remain  bounded  in  the  simula- 
tion output. 

The  equations  for  generating  random  motion  corresponding  to  a 
specified  PSD  will  now  be  developed.  The  specified  PSD,  s(w),  can  be 
generated  by  passing  continuous  white  noise  n(t)  through  a shaping 
filter  H(jm). 

S<u>)  - |»( jto>  |2  sn(w)  (units2/Ha) 
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Figure  4-2.  Random  motion  module  mechanization. 


5 


It  is  convenient  to  make  the  white  noise  dimensionless  and  of  unit 
amplitude.  This  allows  all  of  the  PSD  characteristics  to  be  incorporated 
into  the  transfer  functions  H(j<o).  The  power  spectral  density  for  unity 
dimensionless  white  noise  is 

S (to)  =1  (1/Hz) 

n 

Thus,  the  PSD  of  the  output  of  the  shaping  filter  is 
S (to)  = [H ( jo>)  | 2 (units2/Hz) 

A PSD,  such  as  the  one  in  Figure  4-1,  can  be  approximated  with 
reasonable  accuracy  using  a set  of  second-order  shaping  filters  connected 
in  parallel  as  shown  in  Figure  4-3.  Each  peak  of  the  desired  PSD  is 
generated  by  one  filter.  The  input  white  noise  sources  are  uncorrelated. 


Figure  4-3,  Generation  of  complex  shaped  PSD 


Each  of  the  second-order  shaping  filters  has  the  transfer  function 


Hi  (joi) 


(4-1) 


where 


= undamped  natural  frequency  (rad/s) 
E = damping  ratio 

K = gain  (units) 

The  PSD  of  the  filter  output  is 


The  output  PSD  is  completely  determined  by  K,  <o  , E,  These 

n 

coefficients  can  be  related  to  the  following  characteristics  of  the 
resonant  peak  (see  Figure  4-1) . 

2 

- peak  amplitude  (units  ) 

BW  - half-power  bandwidths  (rad/s) 

- frequency  at  peak  (rad/s) 

The  program  user  specifies  A , BW,  w for  each  of  the  filters. 

P o 

A maximum  of  five  filters  can  be  used  to  approximate  each  PSD. 
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The  program  computes  K,  w , E from  A , BW,  u)  using  the  follow- 

n p o 

ing  relations: 


a = — (0  < a < Jl) 


2 \2*1 1/4 


4 - k - m 


(4-3) 


2 \ 1/2 


ft  - w \ 


(1  - (1  2E2)21A_ 


The  second-order  shaping  filter  in  Eg.  (4-1)  will  be  placed  in 
state  variable  form  to  generate  the  difference  equations  for  computing 
the  sequence  of  numbers  which  simulates  the  random  motion.  The  mechanized 
equations  use  the  following  transform  pair  to  relate  the  PSD  to  the 


auto-correlation  functions : 


+00 

jr(x)q  ^2irfT  dr  (units2/Hz) 


■fco 

R(t)  . - J stfje^21^  df  (units) 


(4-4a) 


where  R(*)  has  the  dimensions  of  (units)  and  S(f)  has  dimensions  of 
2 

(units  /Hz) . The  transform  is  frequently  written  as 
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M/we  :r  c i-il!^jiit?»y<W-f*t<»: 


S(w) 


+00 

/■ 


R(i)e  ^WT  dt  (units2/Hz) 


(4- 4b} 


/ 


R(t)  = / StwJe^*  d in  (units2) 


Notice  that  the  definition  of  the  power  spectrum  S(u>)  has  not  changed 

2 

so  the  dimensions  are  still  (units  /Hz) . Most  PSDs  (such  as  the  PSD  shown 
in  Figure  4-1)  are  plotted  using  the  single-sided  transform 

+00 

SjU)  a 2 


l 


R(x)e 


-jttiT 


dt  + R(0)6(a>) 


m > 0 


That  is,  the  PSD  is  plotted. over  the  positive  frequency  domain  and,  ex- 
cept for  the  dc  component,  is  scaled  by  a factor  of  two.  The  user  should 
scale  the  single-sided  pea)?  amplitude  A^j  the  program  scales  A^  to  the 
two-sided  transform  for  the  internal  calculations.  The  continuous-state 
representation  for  the  shaping  filter  is 

X(t)  * FX(t)  + Cn(t)  (4-5) 

y(t)  «*  HX(t) 

where  n(t)  is  tlie  dimensionless  unit  white-noise  input  to  the  shaping 
filter  and  ¥(t)  is  the  filter  output.  The  scalar  equations  are 


Y(t)  » X^t) 
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(4-6) 


The  mean  and  covariance  of  x(t)  are 


E[X(t)l  = 0 


P(t)  » ^(t]Lt0)  P(tQ)  ♦ (t^)'1  + J <Ht,T)CQCT4>(t,T)' 


(4-7) 


where 


P(tQ) 


» E(X(t)  X(t)  ] 


“ E{X(tQ)  X(tQ)  ] 


E(n(t)n(T)*J  « Q6(t  - t) 


The  magnitude  of  Q follows  from  the  transform  pair. 


S(w)  » / Q6(t.)o  J dt  * Q 


Since  n(t)  is  dimensionless  unity  white  noise 


Q « 1 (l/iiz) 


The  discrete  time  state  representation  for  the  shaping  filter  is 


XCtu...  •■'♦(t..  t ) ♦ BW 


(4-8) 


where  is  a vector  of  dimensionless  discrete  white-noise  samples  of 

unit  magnitude.  The  elements  of  W are  mutually  uncorrelated.  There- 

— n 

fore, 


ElW  WT] 

-tt-m 


Z6 


mn 


The  covariance  equations  for  the  discrete  system  is 

P(tn«)  * ♦‘W.  V Plt„> 

Q = E[BW  WV]  » BBT 
*n  ~n  n 


(4-9) 


The  objective  is  to  simulate  the  continuous  system  (Eq.  (4-5)) 


X(t)  - FX(t)  + <y»(t) 


V(t)  » HX(t) 


with  the  discrete  system  (£q.  (4-8)) 


X *»  $(t  _ , t )X  + BW 

-n+l  r.  n ‘-n 


n 


HX 


such  that  both  systems  have  identical  first  and  second  moments.  The 
elements  of  are  zero  mean  x 
systems  are  zero  mean  because 


elements  of  are  zero  mean  random  numbers  with  unity  variance.  Both 


Eln(t))  - Etwj  - 0 
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To  make  the  covariances  in  Eq.  (4-7)  and  (4-9)  equal,  it  is  necessary 
to  equate 


<J>(t,T)  <J>  (t,T) 


t 

n 


dt 


(4-10) 


Thus  far,  the  elements  of  B are  undefined.  They  will  be  chosen 

so  Eq.  (4-10)  is  satisfied.  The  transition  matrix  over  the  interval 

T = T , , - t is 
n+1  n 


V - ♦(T)  - 


/*11  *12  \ 
'*21  *22 


d)  — 2Eo)  d>  + d) 
911  ' nv12  y22 


1 -EU)„T  . m 
L _ = — e b sin  w_T 

12  w R 

K 


(4-11) 


21 


2. 

•Un*12 


22 


e E(i)nT  cos  a)  T - Eu)  <(>,- 
R n 12 


where 
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Substituting  the  expressions  for  (t , t)  and  t in  Eq.  (4-10)  and 
carrying  out  the  integrations  gives 


Q 


n 


’ll 

’12 


q12 

q22 


- 2Eaj  T r_ 

E n [Eu>  cos  (2o)  T) 
n R 


0)  sin  (2to  T)  ] 


2 2 
4w  (o_ 
n R 


1 - e'2E“nT 


2 2 
4E(o  (o„  4(o  ui 

n R n R 


, -2E(onT  E(0  cos  (2b)  - to„  sin  (2b) 
1 - e " n R 


i22  2 

4EV„„ 


4to, 


i-2E(onT  [Eion  cos  (2(orT  + 2b)  - toR  sin  (2uRT  + 2b)] 

4(o  I 


*12 


-2Ea>nT  , . . cos  b + Euj  sin  b 

3 “'n  - 1 sin  b _R n 


4E(o: 


4(o 


R n 


-*?F  T 

e n [Eto  sin  (2(o_T  + b)  + w_  cos  (2(o_T  + b)  ] 
n R R R 

. 2 
4co„(0 
R n 


where 


b ■ oin  1 (E) 
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The  elements  of  B are  found  by  equating 


2 2 

bn  + bi2 


= Ka)nqil 


2 2 
b12  + b22 


= Ku,nq22 


bllb21  + b12b22 


K“n  qi2 


There  are  three  equations  and  four  b„  coefficients  so  one 
coefficient  is  arbitrary..  Let  b^  = 0,  then 

b22  = V*  “a 


b12  = y/t  <312/Vq22 
bll  ” VK"  “n  Vqn  " qi2/q22 


-13) 


This  gives  all  the  coefficients  needed  for  the  discrete  system 


The  equations  used  to  generate  random  displacement  fix  are 


(4-14) 


where  W^(n)  and  W2(n)  are  numbers  generated  by  a zero  mean,  unity  variance 
Gaussian  random  number  generator.  Random  velocity  is  obtained  by  differ- 
encing consecutive  samples  of  6X. 


6v(n)  * [6X(n)  - 5X(n-l)]/T 


(4-15) 


The  system  is  stable  for  all  damping  ratios  greater  than  zero. 
Thus,  the  covariance  matrix  will  always  approach  steady  state.  Since 
the  random  vehicle  motion  is  assumed  to  exist  prior  to  the  start  of  the 
simulation,  the  initial  conditions  are  chosen  so  the  system  starts  in 
the  stationary  condition.  The  required  relations  can  be  derived  directly 

m 

from  the  continuous  system  covariance  equation. 


P(t)  = FP (t)  + P(t)FT  + CQCT 


In  steady  state,  P(t)  is  zero. 


T T 

0 * FP„„  + P F + Q 

SS  SS 


(4-16) 


The  F and  C matrices  are  given  in  Eq.  (4-6).  Q equals  unity. 


The  steady-state  covariances  are 


Kw 

> = _ n 

11  4E 


12 


(4-17) 


22 


Kco 
n 

4E 


The  initial  state  is  generated  from 


X(t-)  = A W 

0 — o 


(4-18) 


where  consists  of  zero  mean  unity  variance  numbers  from  the  Gaussian 
number  generator  and  the  elements  of  A are  chosen  so  P(t  ) equals  P . 

U So 


E[x(tQ)  2i(t0)  1 


T T 

AE  [W  W ] A 
- o~o 


P„e  = AA 
SS 


c :) 


ail  ai2  \ /ail  a2l\ 

/ \ai2  a22/ 


a21  a22/  \ai2  a22, 


Solving  for  a^_.  terms  gives 


11 


22 


M 

1 4E 


S12  " S21 
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The  initial  state  variables  are  confuted  from 


4.1.2  Random  Angular  Motion 

The  output  of  the  random  angular  motion  generator  is  random  angular 
velocity  in  body  coordinates.  The  average  angular  velocity  over  two 
consecutive  simulation  cycles  is  used  to  compute  the  transformation  from 
nonvibratory  to  vibratory  body  coordinates. 


“ R * 64-l)/2 

c£n  - c£nF  6jj)b 

ft  first-order  attitude  algorithm  is  used  to  update  the  attitude 
matrix  each  computation  cycle 

c£ft(t  + T)  - M(T)c^ 
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The  attitude  matrix  is  orthonormal ized  each  cycle  by  the  following 
equation  which  is  accurate  to  first  order 


c^lortho)  . 4 - ±c£n  [c£n<£n  - l] 
b 

where  C^n  is  the  nonorthonormalization  raatri:;. 

The  total  angular  rate  in  body  coordinates  is  formed  by  adding 
the  random  angular  velocity  <5w  to  the  angular  velocity  obtained  from 
the  flight-path  generator 

of3  «»  tJ3  + 6a!3 


where 

of  » total  angular  velocity  in  body  coordinates 
bn 

» angular  velocity  from  flight-path  generator 
wb  « random  angular  velocity 

The  time  derivative  of  total  angular  rate  is  obtained  by  differ- 
encing consecutive  samples  of  total  angular  rate. 

) ! 

•b  / b b \ j 


4.1.3  Random  Linear  Motion 

The  output  of  the  random  linear  motion  generator  is  random  linear 
velocity  in  body  coordinates  * Random  acceleration  is  computed  by  differ- 
encing consecutive  samples  of  velocity. 


6ab  » 
— n 


(&  - -^-i)/T 


The  total  acceleration  in  the  vibrating  body  frame  is  obtained 

g 

by  adding  the  random  acceleration  6a  to  the  acceleration  obtained  from 
the  flight-path  generator  after  it  has  been  transformed  in  vibrating 
body  coordinates. 


b 

a 


b 

a 


where 

a*  = total  acceleration  in  vibrating  body  frame 
bn 

a = acceleration  from  flight-path  generator 

6a  = random  acceleration 

b 

Ct)n  « transformation  from  nonvibrating  to  vibrating  body  frame 
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SECTION  5 


GYRO  MODELS  AND  COMPENSATION 


The  program  is  structured  so  that  either  SDF  or  laser  gyro  strap- 
down  systems  can  be  simulated.  The  gyro  modules,  gyro  compensation  modules 
and  initialization  data  files  are  designed  to  be  functionally  identical  and 
can  be  replaced  simply  by  specifying  the  correct  data  sets  for  the  desired 
class  of  instruments.  The  dataset  names  are  listed  in  Table  5-1. 


Table  5-1.  Dataset  names  for  SDF  and  laser  gyro  modules. 


DATASET  NAMES 

Source  Programs  and  Data  Files  „„„  _ 

^ SDF  Gyros  Laser  Gyros 


Gyro  model 
Gyro  compensation 


SDFGYRO . FORT  RLGGYRO . FORT 

SDFCOMP . FORT  RLGCOMP . FORT 


Gyro  model  data 
Gyro  compensation  data 


SDFGYRO. DATA 
SDFCOMP. DATA 


RLGGYRO. DATA 
RLGCOMP. DATA 


The  SDF  gyro  and  gyro  compensation  modules  are  described  in  Section  5.1 
The  RLG  laser  gyro  and  gyro  compensation  modules  are  described  in 
Section  5.2. 


i'-’  ttsaM 


5.1  SDF  Gyro  Module 

The  SDF  gyro  module  (gyros)  simulates  three  body  mounted  single- 
degree-of- freedom  gyros  operating  in  torque-to-baiance  modes.  The 
single-degree-of- freedom  gyro  contains  a multitude  of  possible  error 
sources  when  it  is  used  in  strapdown  mechanization  due  to  the  high 
angular  rates  which  may  be  present.  Only  the  most  dominant  of  these 
error  terms — anisoinertia,  cross  coupling,  and  output  axis  accelerations — 
are  included  in  the  model.  Float  suspension  and  float  misalignment  effects 
are  ignored.  The  gyro  model  contains  the  following  error  sources: 

• Anisoinertia 

• Spin-output  axis  cross  coupling 

• Output  axis  acceleration 

• Gyro  input  axis  misalignment 

• Scale-factor  error  (both  positive  and  negative) 

• Scale-factor  rate  sensitivity  (both  positive  and  negative) 

• Gyro  bias 

• Exponentially  correlated  random  drift 

• Turn-on  transient 

• G and  G-squared  coefficionts 

• Angular  quantization 

Tho  program  allows  the  user  to  specify  one  of  three  models  for  the  re- 
balance loop: 

a.  "Performance1*  model  which  ignores  gyro  inertia  and  gyro 
damping. 

b.  A first  order  differential  equation  model  which  contains  gyro 
damping. 

c.  A second  order  differential  equation  model  which  contains 
both  gyro  damping  and  gyro  inertial. 
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The  validity  of  these  various  approximations  will  be  examined  in 
Section  5.1.1. 

Section  5.1.2  derives  the  equations  describing  the  dynamics  of 
the  SDP  gyro.  The  torque-rebalance  loop  is  then  constructed  using  the 
dynamics  of  the  instrument  about  its  output  axis.  Section  5.1.3  de- 
scribes the  mechanization  of  the  SDF  gyro  module. 

5.1.1  Dynamics  of  the  SDF  Gyro 

This  section  develops  the  dynamics  of  the  SDP  gyro  to  the  extent 
they  are  represented  in  the  gyro  module.  While  the  dominant  character- 
istics of  the  gyro  are  included  in  the  model  other  known  mechanisms 
have  been  omitted.  The  following  assumptions  are  implicit  in  the  deri- 
vation t 

a.  The  float  is  perfectly  aligned  with  the  gyro  case  and  can 
rotate  relative  to  the  case  only  about  the  output  axis.  This 
implies  infinite  rotational  suspension  stiffness, 

b.  The  products  of  inertia  of  the  float  assembly  are  zero. 

c.  The  gyro  rotor  gimbal  is  perfectly  rigid  relative  to  the 
float. 

d.  The  gyro  rotor  is  maintained  at  a constant  angular  velocity 
even  in  the  presence  of  angular  vibrations. 

The  development  of  the  gyro  dynamics  will  closely  follow  that  of 
Britting  14] . 

A simplified  illustration  of  the  SDP  gyro  is  presented  in 
Figure  5-1,  where  the  input,  output,  and  spin  (i,  o,  s)  axes  form  an 
orthogonal  set  fixed  to  the  gyro  case.  The  gyro  dynamics  are  obtained 
by  applying  Newton's  second  law  to  the  float  plus  rotor  assembly. 
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Figure  5-1.  Single-degree-of- freedom  gyro. 


where 

“ torque  applied  to  the  float  assembly  the  center 

of  mass 


“ angular  momentum  of  the  float  assembly  about  the 
center  of  mass 


li  j “ ti^e  derivation  of  wrt  inertial  space 


The  time  derivative  of  H_  wrt  the  float  is  more  convenient.  The 

— I 


Coriolis  relation  gives 


" (dt  -f)f  + ^if  * -f 


Expressed  in  float  coordinates  (f) 


f • f / f \ 
~f  - + 


f 


(5-2) 


where  F ( w)  is  the  matrix  form  of  the  vector  product  operator. 


F(uj) 


w. 


“U3  “2 


-u> 


L -0)o  w 0 J 

€*  X 


The  objective  is  to  relate  the  float  dynamics  to  the  angular 
velocity  of  the  gyro  case.  Expressing  the  angular  velocity  of  the 
float  as  the  angular  velocity  of  the  caso  plus  the  angular  velocity 
of  the  float  wrt  the  case  gives 


f it 

(JK,  =*  (1).  + W - 

—if  — ic  — c£ 


_f  , c . c . 
C (w»  ♦ u J 

C “iC  -cf 


(5-3) 


wnere 


-*ic  “ {*V  V Us} 
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According  to  the  initial  assumptions,  the  float  can  rotate  relative 
to  the  case  only  about  the  output  axis;  i.e. 

wC.  = {0,  0 , 0}  (5-4) 

— cf  o 


and  the  principal  axes  of  the  float  assembly  are  perfectly  aligned  with 
the  (i,  o,  j)  case  axes  when  the  float  output  angle  is  zero.  Thus,  the 
transformation  from  case  to  float  is 


(5-6) 


The  angular  velocity  of  the  float  in  float  coordinates  is 


f 

W . „ «3 

“if 


‘ta.  - 0 u>  ' 

1 O s 


w ♦ 0 
o o 


, X O iJ 


(5-7) 


The  angular  momentum  of  the  float  assembly  is  equal  to  the  angular 
Momentum  of  the  float  global  plus  the  angular  momentum  of  tho  rotor. 


f f £ 

IL  " H + K 
— t g fc 


■ Vfe  * “r 


Since  the  products  of  inertia  are  assumed  to  be  zero, 


U 


fm  m 

I,  1 

V - 0 w ” 

• ■ 

0 

i 

i OS 

I 

to  + $ 

+ 

0 

0 

O 0 

T 

w + i w. 

H 

L S J 

L S O iJ 

b ml 

(5-8) 


where 

H “ rotor  3pin  angular  momentum 

^ I , Ig  - principal  float  moments  of  inertia 

Substituting  Eq.  (5-7)  and  (5-S)  into 

mJ  - il  + (5-2) 

gives  the  following  expression  for  the  dynamics  about  the  output  axis. 

% - Vo - hbi  ♦ (is  - ♦ »(it - W ■ “s' 

+ OHu).  + I u (5™9) 

j o o 


The  torque  applied  about  the  output  axis,  Mq,  is  assumed  to 
consist  of 

« 

1,  Torque  from  the  torque  generator,  M 

2,  Viscous  torque  opposing  output  axis  rotation,  C 6 

3,  Gyro  disturbance  torquua, 
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fi  l«atia  fit  I’.'*!1***  “>"■**'" 


That  is, 


M = M - C.0  + M, 

o tg  0 o d 


(5-10) 


Substituting  Eq.  (5-10)  into  Eq.  (5-9)  are  rearranging  terms  gives 


10+ce  3 Hto . + M. 

o o o o i tg 


(ideal  gyro  equations) 


+ (I  - I,)0).(0 
s i l s 


(anisoinertia  torque) 


2 2 

+ Q(I  - I.) (w.  - u,  ) (anisoinertia  coupling) 

S X x 3 


- OHw 


s 


(cross-coupling  torque) 


- I 0) 
o o 


(output  axis  angular  accelera- 
tion) 


+ M. 


(disturbance  torque 


(5-11) 


Rebalance  loops  aro  designed  using  the  ideal  gyro  equations 


I 0 + C 0 » Ha>.  - M, 

oo-oo  x tg 


Sineo  the  other  terms  are  ignored,  they  produce  errors  in  the 
indicated  angular  velocity  and  must  be  compensated  in  the  navigation 
computer.  While  analog  torquing  may  bo  used  in  low  accuracy  applications, 
pulse  torquing  is  currently  used  for  inertial  navigation  systems.  A 
typical  rebalance  loop  is  shown  in  Figure  5-2 j each  output  pulse  repre- 
sents an  increment  of  rotation  about  the  input  axis.  All  of  the  error 
torques  have  boon  lumped  into  for  convenience. 
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TORQUE  GENERATOR 


Figure  5-2.  Pulse-torquing  rebalance  loop. 

Figure  5-3  shows  the  rebalance  loop  mechanized  in  the  gyro  module. 
Notice  that  the  rebalance  torque  is  generated  by  multiplying  the  float 
output  angle  by  the  linear  gains  K where 


K 


K K K 
sg  e tg 


Thus,  the  gyro  module  simulates  an  analog  rebalance  loop  followed  by  a 
delta-modulator  to  generate  the  incremental  angle  output. 


K “ KCrt 

tg  e sg 


Figure  5-3.  Simulated  analog  rebalance  loop. 
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The  rebalance  loop  equations  are  obtained  by  substituting 


£•><  f 


M - -K9 
tg  o 


into  Eq.  (5-11)  and  rearranging  terms. 


10  + C S + [Hu  + (I  - I.)  (tof  - u2.)  + K]0 

oooo  s SiSl  o 


Hu),  - M.  + (I  — X , ) cu . ui 
i d s i a s 


- I 0) 
o o 


Again,  it  is  convenient  to  lump  all  of  the  error  terms  into  the 


distrubance  torque  M^.  This  gives 


10  +C0  +K0  = Ho),  + M 

o o oo  lo  i d 


K.  a Htl)  + (I  - X . ) (tuf  - wj)  + K 

X S S X S X 


(5-12) 


Mj  a M.  + (I  - I.)u>.u>  -1(1) 

d d s iis  co 


The  gyro  module  allows  the  user  to  selectively  solve  for  0q  in 


Eq.  (5-12)  using  one  of  three  models  for  the  rebalance  loop: 


a.  A "performance"  model  which  ignores  gyro  inertia  and  gyro 
damping:  i.e. 


O a (Hw4  + MJ/K, 

O 1 u 1 


b.  a first  order  differential  equations  which  includes  the  gyro 
damping:  i.e. 


Co0o  + Kl9o  - Hai+  Hd 


.-.-Vi 


c.  The  complete  second  order  representation;  i.e 


I.  0 + C 6 + K,6  = Hu).  + Mj 

oo  oo  lo  i d 

Figure  5-4  shows  the  effect  of  approximating  the  second  order  dif- 
ferential with  first  and  zero  order  differential  equations.  Both  approx- 
imations artificially  increase  the  gain  of  the  rebalance  loop  at  high 
frequencies,  but  the  first  order  approximations  is  reasonable — it  has 
negligible  effect  on  the  bandwidth  of  the  loop.  The  zero  order  approxi- 
mation completely  eliminates  the  loop  dynamics  and  thus  produces  an 
infinite  bandwidth. 

5.1.2  SDF  Gyro  Module  Equations 

The  SDF  gyro  module  simulates  the  three  rebalance  loops  sequentially. 
Figure  5-5  shows  the  data  flow  for  one  of  the  platform  axes — the  other  two 
axes  are  identical.  The  equations  indicated  in  Figure  5-5  are  summarized 
in  this  section. 

(1)  Transformation  into  Gyro  Coordinates 

The  acceleration,  angular  rate  and  derivative  of  angular  rate  gen- 
erated by  the  random  motion  module  (ENV)  is  transformed  into  gyro  (i,  o,  s) 
coordinates 


g 

a « 


q B 
Cy  a 

B - 


a B 
CB  & 


•B 

0) 


Ji  *B 
C'  to 

B - 


9 * 92'  93* 
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Gyro  module  computations  for  each  platform  axi 


g , 

where  C0  is  the  transformation  to  (i,  o,  s)  gyro  coordinates.  The  ele- 

g 

ments  of  C’  are  not  programmed  constants— -they  are  elements  in  the  gyro 

B 

initialization  data  set.  The  gyro  input  axis  misalignment  must  be  in- 
corporated into  the  C9  input  data. 

B 


(2)  g and  g-Squared  Drift 

The  drift  generated  by  the  g and  g-squared  coefficients  is  ob- 
tained by  multiplying  each  coefficient  by  the  appropriate  component  of 


D » K.a.  + K a + K a + K.,aT  + K.  a. a 

g ii  oo  ss  iii  ioio 


+ K,  a. a + K aa  + K a 
is  i s os  o s ss  s 


where 

a9  a {a.,  a , a } 

— i o a 

It  n linear  acceleration  sensitivity  (mass  unbalance)  for 

•**  1,1. 

acceleration  along  the  K n gyro  axis 

® , q- squared  sensitivity  (compliance)  for  acceleration  along 
the  j and  k gyro  axes 


(3)  Gyro  Bias 

The  gyro  bias  coefficient  D_  represents  the  constant  disturbance 

B 

torque  applied  to  the  float.  Typically,  there  is  a shift  in  this  coef- 
ficient each  time  the  system  is  turned  on. 
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(4)  Exponentially  Correlated  Random  Drift 


Exponentially  correlated  random  drift  is  used  to  describe  the 
stability  of  the  instrument  during  periods  of  continuous  operation; 
i.e.,  when  there  are  no  power  cycles.  The  amplitude  and  correlation  time 
of  the  random  component  are  obtained  from  power  spectral  density  plots 
of  instrument  test  data.  The  PSD  and  correlation  functions  are  related 
by  the  transform  pair 


^ I 


R(t)  e ■^U)T  dt 


“ f 

2ir  J 


S(w)  e^1  dm 


The  correlation  function  for  exponentially  correlated  random 
drift  is  uniquely  defined  by  the  amplitude  and  break  frequency  of  the 
corresponding  PSD.  Figure  5-6  shows  the  relations  in  tenra  of  a single- 
sided  PSD  plot  which  is  how  most  PSD  data  is  displayed. 


|R(t)  (unlts^) 

L R(t)  = e'|,|/TR 


S(cu)  (unltsVHz) 


a a 

TD  « 1/m 


Figure  5-6.  Autocorrelation  function  and  PSD  for 
exponentially  correlated  drift. 
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The  gyro  module  simulates  an  exponentially  correlated  time  series 

using 

D (n+1)  = D_(n)e“At/TR  + o2  (1  - e"2At/TR)  W_;  D_(0)  = 0 

R R R R 

where 

Dr  ■ sifted  random  drift 
At  = simulation  time  step 
t =>  correlation  time 

R 

2 

o = steady  state  variance  of  the  random  drift 

WR  a zero  mean,  unit  variance  Gaussian  sample  from  the  random 
number  generator 

Since  the  random  drift  is  initiated  to  zero,  the  mean  square 
drift  approaches  the  steady-state  variance  in  an  exponential  manner. 

o2(t)  ■ e2(l  - e~2<'/,TR) 

(5)  Exponential  Turn-On  Transient 

The  spurious  error  torques  generated  at  turn-on  until  the  thermal 
system  stabilizes  is  modeled  as  a simple  exponential 

D « D (0) 

T T 

where 

0^(0)  « initial  amplitude  of  the  transient 

<■  exponential  decay  time  constant 
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(6) 


Anisoinertia  and  OA  Acceleration  Torques 


is 


The  sum  of  the  anisoinertia  and  output  axis  acceleration  torques 


(I  “ I.Jtt.M  - I 0) 
s i i s o o 


where 

I. , I , I = principal  inertias  of  the  float  about  (i,  o,  s) 
i o s 

u>. , W * angular  velocity  of  the  case  about  the  input  and 

i s 

spin  axes 

u ® time  derivative  of  the  angular  rate  of  the  case 

about  the  output  axis 

(?)  Rebalance  Loop  Dynamics 

The  float  output  angle  6q  may  be  computed  using  one  of  the  three 
options  shown  below. 


10  + 0 0^  + K 0 “ M 

o o o o 1 o 


L 


PERFORMANCE  MODEL 

J 


FIRST-ORDER  MODEL 


J 


SECOND-ORDER  MODEL 
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where 


K,  = Hu)  + (I  - I.)  (to2  - w?)  + K 

i s s i s 1 

M = H (to,  - D - D - D_  - D ) + M.. 

i B R T g w 


bias  drift 

exponentially  correlated  random  drift 

turn-on  transient 

g and  g-squared  sensitive  drift 

anisoinsrtia  and  OA  acceleration  torques 

This  is  a direct  mechanization  of  the  rebalance  loop  dynamics  shown  in 
Eq.  (5-12).  The  program  defines  gyro  bias,  random  drift,  g and  g-squared 
drift  with  the  opposite  sign  from  the  normal  conventions.  The  impact  of 
approximating  the  rebalance  loop  with  either  the  “performance”  or  first 
order  models  is  discussed  in  Section  5.1.2.  The  discrete  equations  for 
each  options  are: 


D = 

g 

\ * 


(a)  Performance  Model 

0Q(n+l)  - M(n+1)A1 

(b)  First  Order  Model 

6 (n+1)  « 0 (n)  + I-K.0  (n)  + M(nn)JAt 

O O CIO 


0 (o)  * Hu).  (o)/K 

o 1 


(c)  Second-Order  Model 


0 (n+1)  » Q (n)  + 0 (n)  t 

o o o 

Qn  (n+1)  « 0 (n)  + ~ [-C  0 (n)  - K.0  (n)  + M(n+1)]  T 

o o loo  1 o 

o 

0 (o)  =0 

0 

eo<o)  « H^(o)/K 
(8)  Torquer  Scale  Factor 

The  torquer  scale  factor  consists  of  a linear  term  plus  a term 
proportional  to  input  angular  rate. 

S »■  11  + SO  + si(wi)J 

• 11  ♦ “ + slRI  °o)J 


The  scale  factor  coefficients  may  have  different  values  for  positive 
slid  negative  angular  rates.  The  scale-factor  equations  are: 

s(»  ' [l  + s0<+)  * sl<+)  (fi0o)]  °o-° 

s(->  * t1  * (I  »Jj  V'° 

where 

S0(+j , SO^j  “ positl\f>e  ahd.tiOgat.ive  linear  scale  factors 

Si  » positive  ard  fate  sensitive  scale 

factor*  ; 

• - $-19  . - ; . 


■ -V.- 


The  output  pulses  are  then  converted  into  an  indication  of  the 
angular  rotations  about  the  input  axis  of  the  gyro.  The  output  angle 
is  a direct  indication  of  the  angular  velocity. 


w. 


K -1 
— S 
H 


H qn0 


0 

o 


Thus,  the  measured  rotations  about  the  input  axis  is 


a 

50^  ■ w^A t 

K 

- q*effAt 
A0i  « UQi 

The  navigation  equations  reset  A0i  to  zero  each  time  the  platform  attitude 
matrix  is  computed. 


5.1.3  SPF  Gyro  Compensation 

The  SDF  gyro  componsatior;  modulo  (GCOMP)  computes  the  compensation 
for  the  following  error  sources » 

• Gyro  input  axis  misalignment. 

• Scale  factor  error  (positive  and  negative) . 

• Scale  factor  rate  sensitivity  (positive  and  negative) . 

• Bias  drift. 

• g and  g-squarod  sensitive  drift. 

• anisoinortia  torque. 

• Output  axis  acceleration  torque. 

A flow  diagram  of  the  compensation  module  is  shown  in  Figure  5-7. 
The  compensation  equations  for  each  platform  axis  are  summarised  below. 
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SCALE  FACTOR 
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(1) 


Scale-Factor  Compensation 


The  incremental  platform  rotations  about  each  gyro  input  axi9 
are  compensated  for  linear  and  rate  sensitive  scale  factors.  Separate 
scale  factors  for  positive  and  negative  angular  rates  are  provided. 


S 


1 + SO 


(+) 


+ SI 


(+) 


A8  . > 0 

l — 


S 


1 + SO 


(-) 


+ SI 


(-) 


A9 . < 0 

l 


where 


SO^,  = positive  and  negative  linear  scale  factors 

Sl^+),  Sl^  j = positive  and  negative  rate  sensitive  scale  factors 


(2)  g and  g- squared  Sensitive  Drift  Compensation 

The  compensated  incremental  velocities  from  the  accelerometers 
are  transformed  into  gyro  (i,  o,  s)  coordinates. 


Av‘ 


cg  AvA 

b — 


where 


Avg  =»  {Av, , Av  , Av  } of  Kth  gyro 
— i o s 


Av  » compensated  accelerometer  outputs 


transformation  from  body  to  the  K gyro 


The  drift  produced  by  the  acceleration  sensitive  coefficients 


is 


eu>.  = - (K.  a,  + K a + K a ) 

i 11  o o s s 


2 5 

-{K  a + K.  a. a + K.  a. a + K a a + K a ) 
1!  1 lo  1 O IS  1 s os  os  ss  s 


The  program  computes  corrections  to  the  incremental  rotation 
about  the  gyro  input  axis. 


50  = K Av.  + K Av  + K Av 

g ii  o o s s 


2 2 
+(K. .Av.  + K.  Av. Av  + K.  Av.Av  + K Av  Av  + K Av  ) /At 
11  1 10  i o is  i s os  o s ss  s 


where  K and  K are  the  compensation  coefficients. 

K = linear  acceleration  sensitivity  for  acceleration  along 
th 

the  K gyro  axis 

= acceleration  sensitivity  for  accelerations  along  the 
j and  k gyro  axes 


(3)  Gyro  Bias  Compensation  «. 

Gyro  bias  produces  a constant  error  in  indicated  angular  rate 
along  the  gyro  input  axis. 


E“i  = 'db 

The  bias  compensation  is 

60_  = D At 

B B 
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(4) 


Anisoinertia  Torque  Compensation 

Anisoinertia  torque  produces  the  angular  rate  error 


6u). 

1 


I . )u>.U) 

1 IS 


The  anisoinertia  torque  compensation  is 


60n  ■ - “ (I  - I.)A9.A6  /At 

A H s i l s 

where 

H = rotor  spin  angular  momentum 

I . , I = principal  moments  of  inertia 
is 


(5)  OA  Acceleration  Torque  Compensation 

Output  axis  acceleration  torque  produces  the  angular  rate  error 


(0 

o 


The  OA  acceleration  compensation  is 

I 

fie  ~ foe  (n+l)  - A0  (n)  ] / At 

OA  H o O 


where 

A0  «*  incremental  rotations  about  the  output  axis 
o 

I =*  float  inertia  about  the  output  axis 
o 

H » rotor  spin  angular  momentum 
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(6)  Total  Gyro  Compensation 

The  compensated  rotation  about  the  gyro  input  axis  is 

A0 . '.comp)  = S(A0.)  + (60_  + <50  + 60.,  + 60^,) 

i i B g X OA 

where 

S = scale  factor  compensation 

60_  = bias  compensation 

60^  » g and  g-squared  compensation 

<50  j = anisoinertia  torque  compensation 

^8OA  = acceleration  compensation 

The  three  compensated  gyro  outputs  are  then  transformed  from 
the  gyro  input  axes  to  body  coordinates. 


whore 

AO9  =•  compensated  gyro  outputs  along  the  three  gyro  input  axes 
b 

Cg  “ (nonorthogonal)  transformation  from  the  gyro  input  axis 
to  body  coordinates 

ly 

The  transformation  matrix  is  labeled  QMXS  in  the  program.  The 
elements  of  this  matrix  must  be  consistent  with  the  rotation  matrices 
which  specify  the  (i,  o,  s)  axes  of  each  gyro  relative  to  body  coordinates 
These  rotation  matrices  are  labeled  QGBX,  QGBV,  and  QGBZ.  All  of  the 
matrices  are  obtained  from  the  initialization  dataset  for  the  compensation 
module. 
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5.2  LASER  GYRO  MODULE 


The  laser  gyro  module  (GYROS)  simulates  three  body  mounted  gyros 
plus  their  incremental  angle  detectors.  The  laser  gyro  model  contains 
the  following  error  sources: 

• Gyro  input  axis  misalignment. 

• Scale  factor  error. 

• Scale  factor  turn-on  transient. 

• Gyro  bias  drift. 

• Turn-on  transient  drift. 

• Angular  random  walk. 

• White  angular  noise. 

• Angular  quantization. 

The  model  of  the  laser  gyro  is  discussed  in  Section  5.2.1.  section  5.2.2 
describes  the  mechanization  of  the  laser  gyro  module  (GYROS)  and 
Sectipn  5.2.3  describes  laser  compensation  module  (GCOMP) . 

5*2.1  Principles  of  Operation  and  Model  of  the  Laser  Gyro 

This  section  develops  the  principles  of  operation  and  other 
salient  characteristics  of  the  ring  laser  gyro,  a recent  development 
in  optical  technology  which  shows  promise  of  replacing  many  conventional 
inertial  gyros.  These  principles  lead,  in  turn,  to  the  model  imple- 
mented in  this  module  for  simulating  the  performance  of  these  strap- 
down  laser  gyros. 

Tho  laser  gyro  is  basically  a laser  that  has  three  or  more  re- 
flectors arranged  to  enclose  an  area.  These  throe  mirrors,  in  con- 
junction with  the  light-amplifying  material  in  the  laser  path,  com- 
prise two  optical  oscillators— one  that  has  energy  traveling  clockwise, 
and  one  that  has  energy  traveling  counterclockwise  around  the  same  path. 
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The  optical  length  of  the  path  they  travel  determines  the  frequencies 
at  which  these  oscillations  operate.  When  the  enclosed  ring  is  rotated 
in  inertial  space,  the  clockwise  and  counterclockwise  paths  have  effec- 
tively different  lengths.  The  path  difference  in  these  two  directions 
causes  the  two  oscillators  to  operate  at  different  frequencies,  and 
this  difference  is  proportional  to  the  rate  at  which  the  ring  is  rotating. 
The  readout  of  the  laser  gyro  is  then  accomplished  by  monitoring  the 
frequency  difference  between  the  two  lasers. 

The  foundation  of  laser  gyro  technology  is  the  description  of  the 
phenomenon  by  which  the  path  around  the  ring  can  be  different  for  ob- 
servers (photons)  traveling  with  the  direction  of  rotation  than  for  ob- 
servers traveling  against  the  direction  of  rotation.  In  accordance  with 
the  principles  of  general  relativity,  two  observers,  traveling  around  a 
closed  path  that  is  rotating  in  inertial  space,  will  find  that  their 
clocks  are  not  in  synchronization  when  they  return  to  the  starting  point 
(traveling  once  around  the  path,  but  in  opposite  directions).  The  ob- 
server traveling  in  the  direction  of  rotation  will  experience  a small 
increase,  and  the  observer  traveling  in  the  opposite  direction  will 
experience  a small  decrease  in  clock  time.  The  difference  in  tho  read- 
ings of  these  clocks  depends  upon  the  inertial  rotation  rate,  the  area 
enclosed  by  tho  path,  and  the  speed  of  light.  Since  the  laser  gyro  uses 
photons,  traveling  at  the  speed  of  light,  for  observors,  the  time  dif- 
ference appears  as  an  apparent  length  change  in  the  two  paths  equal  to 
tho  observers*  velocity  times  their  time  difference.  Accordingly,  oven 
though  both  observers  traverso  an  identical  physical  path,  thoy  soo  an 
apparent  length  change  which  is  proportional  to  the  enclosed  area  of 
the  path  and  its  rotation  rate  in  inertial  space. 

Tho  fundamental  boundary  condition  governing  tho  operation  of  tho 
laser  gyro  is  that  the  laser  wavelength  X must  be  equal  to  an  integer 
function  of  the  optical  length  L around  the  cavity.  That  is 

X » h/N  (1) 
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5 6 

where  N is  a large  integer  typically  in  the  range  of  10  to  10  . A 
length  change  of  AL  will  cause  a wavelength  change  of 

AX  = AL/N  (2) 

The  corresponding  frequency  change,  by  logarithmic  differentiation, 
is  given  as 


Af 

f 


(3) 


Accordingly,  given  small  length  differences  AL  and  reasonable 
cavity  lengths  L,  the  operating  frequency  should  be  as  high  as  possible. 

The  relationship  between  inertial  input  rate  Q and  apparent  length 
change  AL  can  easily  be  shown  to  be 


AL 


(4) 


where  A and  C represent  the  enclosed  area  and  the  speed  of  light, 
respectively. 

Thus,  the  equation  relating  ft  to  Af,  in  terms  of  gyro  size  and 
wavelength,  is  determined  by  substituting  Eq.  (4)  into  Eq.  (3),  yielding 


Af 


4Afi 

XL 


(5) 


This  represents  the  ideal  laser  gyro  equation,  but  fails  to 
account  for  the  three  principal  error  sources  effecting  the  operation 
of  all  laser  gyros.  In  descending  order  of  importance,  they  are: 
lock-in,  null  shift,  and  mode  pulling. 
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Lock-in  arises  due  to  the  mutual  coupling  between  the  oppositely 
directed  traveling  waves  and  has  been  the  focus  for  most  of  the  effort 
to  eliminate  the  laser  gyro's  error  sources.  When  the  rotation  rate 
is  reduced  below  some  critical  value  (called  the  lock-in  threshold) , 
the  frequency  difference  between  the  oppositely  traveling  waves  syn- 
chronizes to  a common  value.  Thus,  for  rotation  rates  below  the  lock- 
in  threshold,  the  laser  gyro  is  unresponsive  to  rotations. 

The  principal  method  employed  to  minimize  the  effect  of  lock-in 
is  to  bias  the  laser  gyro  with  an  artificial  input  rotation  such  that 
zero  inertial  rotation  rate  corresponds  to  a point  exterior  to  the  locked 
region.  As  a result,  the  system  operates  at  all  times  with  a differential 
path  and  a frequency  difference  with  respect  to  the  oppositely  traveling 
beams.  Such  a bias  can  be  introduced  by  any  technique  which  introduces 
an  anisotropy  in  the  index  of  refraction,  again  with  respect  to  the 
oppositely  directed  beams,  and  is  commonly  effected  by  either  mechanical 
or  optical  means,  in  turn,  each  of  these  general  techniques  can  be  em- 
ployed in  two  ways:  by  holding  the  bias  fixed,  measuring  the  input  rate 
or  angle,  and  observing  differences  in  the  laser  gyro  output  rate  due 
to  the  input  rate,  or  by  introducing  a negative-to-positive  alternating 
bias. 

Due  to  a reduction  in  the  requirements  on  the  absolute  stability 
of  the  magnitude  of  the  bias,  alternating  bias  techniques  are  generally 
employed.  The  idea  behind  this  method  is  to  minimize  the  (switching) 
time  that  the  laser  spends  in  the  locked  region.  Since  the  laser  gyro 
is  basically  an  integrating  rate  instrument,  only  the  not  rotation  ap- 
pears  in  the  output.  The  basic  premise  in  the  alternating  bias  technique 
is  that  the  bias  is  Perfectly  symmetrical.  Thus,  the  laser  is  allowed 
to  remain  unlocked  over  most  of  the  "dither"  cycle,  and  hence  the  gyro 
is  responsive  to  any  nonzero  input  rate. 
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The  error  generated  when  the  laser  gyro  is  biased  by  a sinu- 
soidal dither,  at  a rate  several  orders  of  magnitude  greater  than  the 
lock-in  rate,  is  characterized  by  a one-dimensional,  random-walk  pro- 
cess. Although  the  mean  or  expected  error  is  zero,  the  root  mean  squ< 
error  is  not  and  is  given  by 

rms  error  (t) 

where 

K =>  gyro  scale  factor  (arcsec/count) 

t a time  (s) 

a lock-in  threshold  (deg A • arcsec/s) 

^ « peak  dither  rate  (deg/h  » arcsec/s) 


Numerical  values  are  a function  of  the  gyro  size  and  operative  parameters. 
As  an  example , the  Honeywell  GG/1300  laser  gyro  (5.7  inch,  visible)  has 
K *•  1.5  arcsec/count, 

« 3.77  x 105  deg/h 

150  Hz) . Accordingly,  for  this  particular  instrument 

1/2 

rms  error  (t)  - 1.01  t arcsec 


0^  typically  0.5  deg/s  (1800  deg/h),  and 
(amplitude  of  400  arcseconds  at  a frequency  of 


In  addition  to  any  errors  inherent  in  the  laser  gyro,  there  is  additional 
uncertainty  associated  with  the  sampling  process. 

With  regard  to  other  error  sources,  a null  shift  (bias)  is  intro- 
duced into  the  output  of  the  laser  gyro  by  any  effect— '"exclusive  of  iner- 
tial rotation  rates— which  makes  the  clockwise  optical  path  length  differ- 
ent from  the  counterclockwise  path  length.  These  effects  arise  from  the 
anisotropy  (nonreciprocation)  of  the  optical  parameters  within  the  cavity 
with  respect  to  radiation  traveling  in  two  directions  and  cause  a shift 
of  the  beat  frequency  (output)  versus  rotation  rate  ideal  curve.  Unless 
the  gyro  is  properly  designed,  null-shift  errors  can  be  orders  of  magni- 
tude greater  than  input  rates. 
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On  the  other  hand,  mode  pulling,  more  correctly  known  as  anomalous 
dispersion  correction,  is  a phenomenon  which  causes  variations  in  the  gyro 
scale  factor.  The  difficulty  in  maintaining  scale  factor  accuracy  does 
not  lie  in  the  geometric  dimensions,  but  rather  in  the  index  of  refraction. 
The  anomalous  dispersion  due  to  the  atomic  transitions  in  the  lasing 
medium  creates  a gain-dependent  correction  to  the  scale  factor  which  is  a 
strong  function  of  the  operating  characteristics  of  the  laser. 

5,2.2  Laser  Gyro  Model  Equations 

The  laser  gyro  module  (GYROS)  simulates  three  strapped-down  laser 
gyros.  The  gyro  module  computations  are  shown  in  Figure  5-8.  The  equa- 
tions are  summarized  in  this  section.  The  equations  for  each  platform 
axis  are  identical. 

(1)  Scale-Factor  Turn-On  Transient 

The  gyro  scale  factor  consists  of  a bias  plus  an  initial  turn-on 
transient.  The  turn-on  transient  is  modeled  as  a simple  exponential 
function. 

\ « ST(0)e“t/TS 

where 

5^(0)  “ initial  scale  factor  transient  amplitude 

t_  <•  exponential  transient  time  constant 
& . 

(2)  Scale -Factor  and  Gyro  I h Misalignment 

The  throe  gyro  input  axes  are  aligned  with  the  body  axes  except 
for  the  gyro  misalignment  angles  which  are  relatively  small.  The 
transformation  from  body  to  gyro  I A axes  is 

CJ?  * X + 5C? 

B o 

The  program  combines  the  gyro  input  axis  misalignments  with  the  total 
scale  factor  errors. 
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SCALE  FACTOR 
TURN-ON 
TRANSIENT 
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V21 

'31 

K32 

V 32 

"V31 

V 

The  elements  on  the  main  diagonal  are  total  scale  factor  error. 
S 


3 + S 

0 T 


SQ  « constant  scale  factor  error 


S^.  - exponential  turn-on  transient 


The  off-diagonal  elements  are  the  input  axis  misalignments. 
ViD 

platform  axis 


» misalignment  of  the  itli  gyro  input  axis  about  the 


The  laser  gyro  modulo  obtains  the  values  of  the  elements  from 
the  module  initialisation  file.  The  elements  are  specified  directly 
rather  than  the  values  for  scale  factor  and  1A  misalignments.  The  ele- 
ments on  the  main  diagonal  are  the  values  of  constant  scale  factor  error. 


K,  . - S 

ii 


<pp«) 


The  off-diagonal  elements  correspond  to  the  input  axis  misalign- 
ments except  for  sign. 


K,  , **  sv_  {microrad) 

1 J 


Thu  program  then  adds  the  scale-factor  transients  to  the  main 
diagonal  elements  to  form  the  complete  rate  sensitive  matrix. 
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(3)  Gyro  Bias 


The  gyro  bias  D0  produces  a constant  error  in  the  indicated  angu- 
lar velocity  along  the  gyro's  input  axis.  Typically  there  is  a shift 
in  this  coefficient  each  time  the  system  is  turned  on. 

(4)  Exponential  Drift  Turn-On  Transient 

The  turn-on  drift  transient  is  modeled  as  a simple  exponential. 

Dm  = Dm(0)e”t/TD 
T T 


where 

0^(0)  = initial  amplitude  of  the  drift  transient 

td  = exponential  transient  time  constant 


(5)  White  Angular  Noise 

The  measured  angular  rotation  about  the  gyro  input  axis  is  corrupted 
by  angular  random  walk  and  angular  noise  which  is  essentially  uncorrelated 
with  time.  That  is,  the  correlation  function  of  the  noise  is 


VT> 


where  t is  of  the  order  of  several  milliseconds  or  less.  The  program 
N 

approximates  the  noise  process  as  discrete  white  noise  since  the  correla- 
tion time  is  substantially  smaller  than  the  simulation  time  step.  The 
discrete  white  noise  model  is 

J(VntU  " °nwr 


where 

2 

a - variance  of  the  noise 
N 

W - zero  mean,  unit  variance  Gaussian  sample  from  the  random 
n 

number  generator 


(6)  Angular  Random  Walk 


The  amplitude  of  a random  walk  process  is  obtained  from  its  power 
spectral  density  or  by  measuring  its  mean  square  growth  rate  as  illustrated 
in  Figure  5.9. 

The  random  walk  process  is  simulated  using  the  difference  equations 


5<V"+1>  - SeRW(n>  + <‘44t)1/2  wn 

where 

2 

a = amplitude  of  the  random  walk  process 
RW 

At  = simulation  time  step 

W = zero  mean,  unit  variance  sample  from  the  random  number 

n 

generator 


Figure  5-9.  Mean  square  growth  rate  and  PSD  for  random  walk 
process. 
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(7)  Measured  Rotation  Angle 


The  total  measured  rotation  angle  0_^  about  the  gyro  input  axis  is 


a) 


B 


(I  + K.,)  OJ  + D + D„ 

w —ib  b a 


e.  = mAt  + 66rn , + 66„, 
—i  — WN  RW 


where 


co 


B 


■IB 


= angular  rate  of  the  body  in  body  coordinates 


= gyre  scale  factor  and  IA  misalignment  matrix 


D 


B 


= gyro  bias 


T 


WN 


= turn-on  drift  transient 


= angular  white  noise 


60 


RW 


= angular  random  walk 


At  = simulation  time  step 


(8)  Quantized  Output  Angle 

The  pulse  generator  computes  the  number  of  output  pulses  produced 
by  the  continuous  angle  0^  plus  the  stored  angle  residual.  The  output 
pulses  are  then  converted  into  a quantized  indication  of  angular  rotation 
A0^  about  the  input  axis  of  the  gyro. 


0 

= 0. (n+1)  + R(n) 

l 

n0 

= Integer  ( 0/q) 

60. 

l 

= qnQ 

R.(n+1) 

“ 0 - q°e 

A0i(n+1) 

■ A0. (n)  + qn 

i 0 
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where 


0^  = continuous  indication  of  rotation  about  gyro  input  axis 

R = quantized  angle  residual 

n^  = number  of  pulses 

A8^  = quantized  indication  of  rotation  about  input  axis 
q = angular  quantization 

The  navigation  equations  reset  A8^  to  zero  each  time  the  platform  attitude 
matrix  is  computed. 

5.2.3  Laser  Gyro  Compensation 

The  laser  gyro  compensation  module  (flCOMP)  accepts  incremental 
A0  rotations  from  the  gyros  and  produces  compensated  A0  rotations  in 
body  coordinates.  The  compensation  module  computes  the  compensation 
for  the  following  error  sources: 

• Gyro  input— axis  misalignment 

• Scale- factor  error 

• Scale-factor  turn-on  transient 

• Gyro-bias  drift 

• Turn-on  transient  drift 

The  indicated  angular  rate  when  scale-factor  errors,  IA  misalignments, 
and  gyro  drifts  are  present  is 

4 - <1  + V“ib  + 2b  + 2t 
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where 


= gyro  scale  factor  and  IA  misalignment  matrix 
Kg  = gyro  bias 
Dt  = turn-on  transient  drift 

Hie  estimate  of  platform  angular  rate  given  and  estimates  of  the 
error  parameters  are 

%b  * 11  + V"1  <2i  --2b  -£t> 

* « - V (-l  - .5b  ■ -V  * 0 

Since  the  output  of  the  strapdown  gyrus  is  incremental  rotation 
angle  instead  of  angular  rate,  the  compensation  is 

60  = A0.  - (D  + D) At 
— -i  '-B  — T 

A®0  " ” (1^60) 

where 

A0  » compensated  incremental  rotation  in  body  coordinates, 
c 

A0i  = incremental  rotation  obtained  from  the  gyros. 

D_  = gyro  bias  compensation. 

5 

Dt  » turn-on  transient  drift  compensation. 

Ky  - scale-factor  and  IA  misalignment  compensation. 

At  = platform  attitude  matrix  update  interval. 


The  off-diagonal  elements  contain  the  compensation  for  gyro  input- 
axis  misalignment.  As  in  the  laser  gyro  simulation,  the  elements  of 
are  specified  directly  in  the  compensation  module  input  data  file 
rather  than  specifying  the  scale-factor  values  and  the  misalignment 
angles 


K.  . 
XI 

" soi 

K.  . 
ID 

~Uik 

(ppm) 

(microrad) 
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(1)  Exponential  Turn-on  Drift  Transient 


A simple  exponential  is  used  to  compensate  the  turn-on  drift 
transient. 

D = D (0)  d~t/,T0 
T T 

where 

DT(0)  = compensation  amplitude. 

= exponential  transient  time  constant. 


(2)  Scale-Factor  Turn-On  Transient 

The  compensation  for  the  scale-factor  transient  is  the  exponen- 
tial function 


. _ -t/x 

Sm  a S.n.  e s 
T T{0) 


where 


T(0) 


- initial  compensation  aplitude. 

--  exponential  transient  time  constant. 


The  elements  on  the  main  diagonal  of  the  matrix  contain  the 
total  scale- factor  compensation 


K. , = SA  + Sm 

xi  0 T 
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SECTION  6 


ACCELEROMETER  MODELS  AND  COMPENSATION 


The  program  simulates  three  body  mounted  SDF  pendulous  acceler- 
ometers including  the  change  in  specific  force  between  the  accelerometer 
locations  and  a specified  location  in  the  vehicle.  The  accelerometer 
compensation  module  contains  compensation  for  the  lever-arm  effects. 

6.1  Accelerometer  Module 

The  accelerometer  module  (ACCEL)  simulates  three  strapped-down, 
SDF  , pendulous,  floated  accelerometers  plus  accelerometer  lever-arm 
effects.  The  output  of  the  module  is  incremental  velocity  in  accelero- 
meter coordinates.  The  SDF  pendulous  accelerometer  contains  numex'ous 
error  sources  in  strapdown  environments  due  to  the  high  angular  rates 
which  are  present.  Only  the  most  dominant  of  these  error  terms  - aniso- 
inertia  and  output  - axis  accelerations  - are  included  in  the  accelero- 
meter model.  Float  suspension,  float  misalignment,  and  products  of 
inertia  are  ignored.  The  SDF  accelerometer  model  contains  the  following 
error  sources: 

• lever-arm  effects 

• anisoinertia 

• output-axis  acceleration 

• accelerometer  input-axis  misalignment 

• scale-factor  error  (both  positive  and  negative) 

• scale-factor  rate  sensitivity  (both  positive  and  negative) 

• accelerometer  bias 
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• exponentially  con elated  accelerometer  noise 

• g and  g-squared  coefficients 

• quantization 

The  program  allows  tha  user  to  specify  one  of  three  models  for 
the  accelerometer  rebalance  loops 

(a)  A "performance"  model  which  ignores  accelerometer  float 
inertia  and  accelerometer  damping. 

(b)  A first-order  differential  equation  model  which  contains 
the  accelerometer  damping „ 

(c)  A second-order  differential  equation  model  which  contains 
both  accelerometer  damping  and  inertia. 

The  validity  of  these  approximations  was  examined  previously  for 
the  SPF  gyro  in  Section  5.1.2.  The  zero-order  and  first-order  approxi- 
mation artifically  increase  the  gain  of  the  accelerometer  rebalance  loop 
at  high  frequencies.  The  first-order  approximation  is  reasonable  because 
it  has  negligible  effect  on  the  bandwidth  of  the  loop  (see  Figure  5-4) . 
The  zero-order  "performance"  model  completely  eliminates  the  loop 
dynamics  , and  thds  produces  an  infinite  bandwidth. 

Section  6.1.1  derives  the  equations  describing  the  dynamics  of 
the  SDF  pendulous  accelerometer.  The  torque-rebalance  loop  is  then 
constructed  using  the  dynamics  of  the  instrument  about  its  output  axis. 
Section  6.1.2  describes  the  mechanization  of  the  accelerometer  module. 

The  accelerometer  compensation  is  discussed  in  Section  6.1,3. 

6.1.1  Dynamics  of  the  SPF  Pendulous  Accelerometer 

The  SDF  floated  pendrlous  accelerometer  contains  a pendulous  mass 
housed  in  a gimbal  which  is  floated  and  damped  in  a viscous  neutral 
buoyancy  fluid.  The  unbalanced  float  makes  the  instrument  sensitive  to 
linear  accelerations.  The  motion  of  the  float  about  the  output  axis  is 
a direct  indication  of  specific  force  along  the  input  axis  of  the 
instrument.  Thus,  the  instrument  may  be  visualized  as  a SDF  gyro  with 
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the  rotor  replaced  with  a pendulous  mass.  Figure  6-1  is  a simplified 
sketch  of  the  device.  The  instrument  coordinates  are  called  the  input, 
output  and  pendulous  (i,  o,  p)  axes. 

/ 


PENDULOUS 
AXIS  (p)  • 


TORQUE 

GENERATOR 


Figure  6-1.  SDF  pendulous  accelerometer. 
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This  section  develops  the  dynamics  of  the  SDF  pendulous  acceler- 
ometer to  the  extent  they  are  represented  in  the  accelerometer  module. 
While  the  dominant  characteristics  of  the  instrument  are  included  in  the 
model # other  known  mechanisms  have  been  emitted.  The  following  assump- 
tions are  implicit  in  the  derivations 

(a)  The  float  is  perfectly  aligned  with  the  accelerometer 
case  and  can  rotate  relative  to  the  case  only  about  the 
output  axis.  This  implies  infinite  rotational  suspension 
stiffness, 

(b)  The  products  of  inertia  of  the  float  assembly  (gimbal  with 
pendulous  unbalance)  are  zero. 

The  development  of  the  accelerometer  dynamics  will  closely  follow 
the  development  in  Section  5.1.1  for  the  single-degree-of-freedora  gyro. 
This  approach  allows  the  two  types  of  instruments  to  be  readily  compared. 
The  dynamics  of  the  accelerometer  are  obtained  by  applying  Newton's 
second  law  to  the  float  assembly.  While  gyros  are  designed  to  be 
perfectly  balanced  about  the  output  axis#  this  is  not  true  for  the 
pendulous  accelerometer.  A known  unbalance  is  deliberately  designed 
into  the  float.  The  accelerometer  dynamics  are  obtained  by  taking 
moments  about  point  (0)  in  Figure  6-1.  Point  (0)  is  fixed  in  the  float 
and  lies  at  the  origin  cf  the  float  (i#  o#  p)  coordinate  system. 

^fo  “ (dt  ^fo)2  * m(dt  — o 14  (dt-jj  (6-1) 

torcue  applied  to  the  float  assembly  about  point  (0). 

angular  momentum  of  the  float  assembly  about  point  (0). 

mass  of  the  float  assembly, 
position  of  point  (0). 

position  vector  from  point  (0)  to  the  center  of  mass, 
time  derivative  of  tt£o  wrt  inertial  space. 


whore 
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Since  point  (o)  is  fixed  in  the  float,  the  angular  momentum 
about  point  (o)  is 


H.  « I 
— fo  o 


0)  -l±r  ) 

-if  \dt  / 


x mr 


(6-2) 


Substituting  Eg.  (6-2)  into  Eg.  (6-1)  gives 


%>  ’ {k  Hfo)  ♦ »(&  i»)  * (<£  4 


_d_ 

dt 


2 

(l  w \ - (-3—  r j x mr 

\ °”lf/i  \dt2  VI  ” 


The  torque  applied  to  the  float  consists  of  mass  attraction  torque  plus 
other  mechanisms. 


i^o 


Sft>  + 


mr  x <jm 


Thus,  tho  reaction  to  the  nongravitational  torques  is 


(6-3) 


where  a is  the  st>eci£ic  forco.  It  is  more  convenient  to  work  with  the 
time  derivative  with  respect  to  the  float.  Hie  Coriolis  relation  gives 
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The  inertia  matrix  is  constant  because  point  (o)  is  fixed  in  the 


float  assembly.  Expressed  in  float  coordinates  (f) , 


= xcAf + p(4f)(io4f)  + “p(^f)^f 


(6-4) 


where  P(.)  is  the  matrix  representation  of  the  vector  product  operator. 


P(w) 


-UL  <0 

3 2 


“3  0 

■"i  “i  0 


The  objective  is  to  express  the  float  dynamics  in  terms  of  the 
angular  velocity  and  specific  force  in  case  coordinates.  Expressing 
the  angular  velocity  of  the  float  as  the  angular  velocity  of  the  case 
plus  the  angular  velocity  of  the  float  with  respect  to  the  case  gives 


J,  - a?  *»' 

-af  -ic  -of 


. c . c , 
C (w,  ♦ a!  A 

e -ic  -of 


(6-S) 


where 


c 

tu. 

—ic 


<“i<  “o'  V 


According  to  the  initial  assumptions,  the  float  can  rotate  rela- 
tive to  the  case  only  about  the  output  axis;  i.e. , 


wCc  * {o,  0 , o) 

— cf  o 


(6-6) 
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and  the  principal  axes  of  the  float  assembly  are  perfectly  aligned  with 
the  (i,  o,  p)  case  axes  when  the  float  output  angle  is  zero.  Thus,  the 
transformation  from  case  to  float  is 
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0 10 


0 *J 


The  angular  velocity  of  the  float  in  float  coordinates  is 


■u>.  - 0 to  1 

i op 


Wj,  » (D  + 0 

—if  ! o o 


to  + 0 to.  I 

> P • OiJ 


The  specific  force  in  float  coordinates  4s 


ra.  - 6 a 
i op 

f f c 

a - C a * a ■ 

— o — * o 


a + 0 a, 

. p or 


Substituting  Eq.  (6-8)  and  Eg.  (6-9)  into 


4.  • v|f  * rl4t’  (1o &>  * 


(6-7) 


(6-8) 


(6-9) 


(6-4) 


gives  the  following  expression  for  tits  dynamics  about  the  output  axis. 
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M = 1*6*  + (I.  - I )M .u  + 6 (I.  - I )((D?  - uj) 
o oo  i pip  oi  pi  P 


+ I w - mr  a.  + 9 mr  a 
o o pi  o p p 
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The  output  axis  equation  was  derived  assuming  the  products  of 
inertia  are  zero;  i.e. , 


I 


!\ 0 0 


0 I 0 

o 

0 0 I 


and  the  pendulous  mass  lies  along  the  negative  (p)  axis 


£ « {o,  o,  -r  } 

***  P 


The  product  rur 


P 


is  called  the  pendulosity. 
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mr 


The  torque  applied  about  the  output  axis,  Mq,  is  assumed  to 
consist  of 

(a)  Torque  from  the  torque  generator,  Mfcg. 

* 

(b)  Viscous  torque  opposing  the  output  axis  rotation,  C^O^. 

(c)  Gyro  disturbance  torques,  K^, 

That  is. 


H 


M 


tg 


C } 
o o 


+ M, 
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Substituting  Eq.  (6-11)  into  Eq.  (6-10)  and  rearranging  terms  gives 


I 8_  + C 8 
o o o o 


Pa,  + M. 
i tg 


(ideal  accelerometer) 


+ (I  - I.) co  u (anisoinertia  torque) 

P r * p 

2 2 

+ (anisoinertia  coupling) 

# 

- IqWq  (output-axis  angular  acceleration) 

- (cross* coupling  torque) 

+ (disturbance  torques) 

(6-12) 

The  accelerometer  rebalance  loop  Is  designed  using  the  ideal  acceler- 
ometer equation 


• * * 
10  + C U 

o o o o 


Pa,  + M. 
i tg 


Since  the  other  terms  are  ignored,  they  produce  errors  in  the 
indicated  specific  force  and  must  be  compel. sated  in  the  navigation  com- 
puter, Current  inertial  navigation  systtons  use  digital  accelerometer 
rebalance  loops.  A typical  digital  loop  is  shown  in  Figure  6-2.  Each 
output  pulse  represents  an  increment  of  velocity  (integral  of  specific 
force)  along  the  input  axis.  All  the  error  torques  have  been  lumped 
into  for  convenience. 

The  rebalance  loop  mechanized  in  the  accelerometer  module  is 
shown  in  Figure  6-3.  The  rebalance  torque  is  generated  by  multiplying 
the  float  output  angle  by  the  linear  gain  K where 


K = 


Ksg  Ktg 


thus,  the  program  simulates  an  analog  rebalance  loop  followed  by  a 
delta-modulator  to  generate  the  incremental  velocity  output. 


Figure  6-2.  Pulse-torquing  accelerometer  rebalance  loop. 


Figure  6-3.  Simulated  analog  accelerometer  rebalance  loop. 
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The  rebalance  loop  equations  are  obtained  by  substituting 


-KG 

o 


into  Eq.  (6-12)  and  rearranging  terms. 


I V + C 5 + [Pa  + (I  - I.)  (“  - to.) 

oo  oo  1 p p i p i 


+ K]0  = Pa.  + M.  + (I  - I. ) (O.to  - id) 

o i d p i ip  oo 


Again,  it  is  convenient  to  lump  all  of  the  sensor  errors  into  the 
disturbance  torque  This  gives 


10  +C0  +K.0 
o o o o 1 o 


c Pa.  + Mi 
i d 


® Pa  + (I  - I.)  <w 
P prp 


-•*> 


+ K 


Md 


M . + (I  ~ I,  ) « - I w 

d p i i p o o 


(6-13) 


The  accelerometer  modulo  allows  the  user  to  selectively  solve  for 
0Q  in  Eq.  (6-13)  using  one  of  three  models  for  the  rebalance  loops 

(a)  A ‘'performance”  model  which  ignores  float  inertia 
and  damping;  i.e.; 

6o  - (Pai  * Md>/Kl 


(b)  A first-order  differential  equation  which  includes  the 
accelerometer  damping;  i.e., 


CO.  K,0 
O 0 1 o 


Pa.  + Mi 
i a 
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(c)  The  conplete  second-order  model?  i.e.. 


I 8 + C 0 

o o o o 


+ K. 


1 o 


Pa,  + M' 
i d 


The  effect  of  using  the  zero-order  or  first-order  approximations 
has  been  discussed,  but  to  reiterate  - while  the  first-order  model  is 
reasonable  the  "performance"  model  produces  an  infinite  rebalance  loop 
bandwidth. 


6.1.2  Accelerometer  Module  Equations 

The  accelerometer  module  simulates  the  three  accelerometer 
rebalance  loops  sequentially.  Figure  6-4  shows  the  computations  for 
one  platform  axis.  The  computations  for  each  of  the  axes  are  identical. 
The  equations  indicated  in  Figure  6-4  are  summarized  in  this  section. 

(1)  Lever-Arm  Effects 

The  displacement  of  each  accelerometer  from  the  vehicle  reference 
location  is  specified  in  the  accelerometer  module  input  data  file.  The 
lever  arm  is  specified  in  body  coordinates.  The  displacement  of  the 
accelerometer  from  the  reference  location  changes  the  output  of  the 
instrument.  The  change  in  specific  force  can  be  derived  from 
Figure  6-5. 

The  specific  force  at  the  accelerometer  location  expressed  in 
body  coordinates  is 
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Figure  6-5,  Geometry  of  the  lever=>arm  effect. 


The  first  and  second  derivatives  of  r are  zero  because  the  lever 
am  is  assumed  to  be  fixed  in  the  body.  Thus,  for  a rigid  lever  arm, 
the  specific  force  at  the  accelerometer  location  is  equal  to  the  specific 
force  at  the  reforonce  location  plus  the  eentropital  and  tangential 
acceleration  produced  by  vehicle  angular  motion. 


ab 

-a 


b 

2b 


, b 

(%b 


rb> 


. *b  b 


(6-14) 


The  program  mechanizes  Eq.  (6-14). 

(2)  Transformation  into  Accelerometer  Coordinates 

The  specific  force#  angular  velocity,  and  derivative  of  angular 
velocity  generated  by  the  random  motion  module  (ENV)  are  transformed 
into  accelorometer  (i#  o,  p)  coordinates. 
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a 


i 


b 

a 


A.  , 
C 1 “ 

D — 


(i  = 1,  2,  3) 


A. 

• l 


03 


•b 

to 


i th 

where  C,  is  the  transformation  to  (i,  o,  p)  coordinates  of  the  i 

° Ai 

accelerometer.  The  elements  of  are  not  programmed  constants — they 
are  specified  by  the  accelerometer  initialization  data  set.  The 
accelerometer  input-axis  misalignment  must  be  incorporated  into  the 
input  data. 


(3)  g and  g-Squared  Coefficients 

The  accelerometer  error  generated  by  the  g and  g-squared  coeffici 

ents  is  obtained  by  multiplying  each  coefficient  by  the  appropriate 
. \ 

component  of  a 1 


a 

”9 


K a Ha  + K,,a.  + K,  a.  a 
o o p p ii  1 io  i o 


K . a , a + K a a + K a" 
ip  i p op  o p pp  p 


where 

A. 


*k 


es 


linear  acceleration  sensitivity  for  acceleration 
til 

along  the  k accelerometer  axis. 


g-squared  sensitivity  (compliance)  for  acceleration 
along  the  j and  k accelerometer  axes. 


6-15 


(4)  Accelerometer  Bias 


The  accelerometer  bias  B represents  the  constant  disturbance 
torque  applied  to  the  accelerometer  float.  Typically,  there  is  some 
shift  in  this  coefficient  each  time  the  system  is  turned  on. 

(5)  Exponentially  Correlated  Random  Noise 

Exponentially  correlated  random  accelerometer  noise  is  used  to 
describe  the  stability  of  the  following  tum-on.  The  amplitude  and 
correlation  time  are  obtained  from  power  spectral  density  plots  of 
instrument  test  data.  The  two  parameters  are  uniquely  defined  by  the 
amplitude  and  break  frequency  of  the  PSD.  Figure  6-6  shows  the  required 
relations.  A more  detailed  discussion  of  an  exponentially  correlated 
random  process  is  presented  in  Section  5.1. 3 (4). 


R(t)  (UNITS2)  S(w)  (UNITS2/Hz) 


Figure  G-6.  Autocorrelation  function  and  PSD  for 
exponentially  correlated  random  noise. 
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The  accelerometer  module  simulates  an  exponentially  correlated 
time  series  using 


-At/r  rn  -2At/t  \ 

aR{n  + 1)  = aR(n)e  + J a 1 1 - e )wn»  aR(o)  * 0 

where 

a simulated  random  accelerometer  noise. 

K 

At  =>  simulation  time  step. 
tr  » correlation  time. 

2 

a & steady-state  variance  of  the  random  process. 

a zero-mean,  unit-variance  gaussian  sample  from  the 
random-number  generator. 


Since  the  random  process  is  initialized  to  zero,  the  mean  square 
value  approaches  the  s toady-state  variance  in  an  exponential  maimer. 


(6)  Anisoinertia  and  Oh  Acceleration  Torques 

The  sum  of  the  anisoinertia  and  output-axis  acceleration  torques 
is 


where 


1, 

i 


) 


(a  u) 


I u 
o o 


principal  inertias  of  the  float  about  (i,  o,  p). 


co 


p 


U 


angular  velocity  of  the  case  about  the  input 
and  spin  axes. 


« time  derivative  of  the  angular  Velocity  of  the 
case  about  the  output  axis. 
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(7)  Rebalance  Loop  Dynamics 


The  float  output  angle  6q  may  be  computed  using  one  of  the  three 
options  shown  below. 


10  + C 0 + 

o o o o 


K. 0 
1 o 


M 


I I 

PERFORMANCE  MODEL 

I I 

FIRST-ORDER  MODEL 

L — j 

SECOND-ORDER  MODEL 

Kt  a Pa  + ( I - I.Ku)2  - w2)  + K 

1 p p i p i 

M “ P(a.  - B - aB  - a ) + H , 
i U q w 


where 

D « accelerometer  bias. 

af,  “ exponentially  correlated  accelerometer  noiso. 

a - g and  g-squared  sensitive  error. 

£ 

Mw  **  anisoinertia  and  output-axis  acceleration  torques. 

This  is  a direct  mech&nizati  h of  tho  rebalance  loop  dynamics 
shown  in  Eq.  (6-13) . The  impact  of  approximating  the  rebalance  loop 
with  either  the  "performance*  or  first-order  models  is  discussed  in 
Section  6.1.  The  discrete  equations  for  each  option  are: 
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(a) 


Performance  Model 


6 (n  + 1)  = M(n  + 1)/K, 

o 1 

(b)  First-Order  Model 

eQ(n  + 1}  = 0Q(n)  + ~-t-K10o(n)  + M(n  4-  l)]At 

o 

0 (oj  = Pa. (o)/K 
o i 

(c)  Second-Order  Model 


6 (n  + 1)  = 0 {n)  + 0 (n)  T 

o o o 

6 {n  4-  1)  » 0 (n)  4-  [— C 0 (n)  - K 8 (n)  + M{n  4-  l))At 

o o loo  1 o 

o 

0 (o)  »0 

o 

0Q(o)  » Pai(o)/K 


The  program  defines  accelerometer  bias,  random  a ecu  le  rosso  to  r~ 
noise,  and  the  v and  g-aqu arad  errors  with  the  opposite  sign  from  the 
normal  convention, 

(8)  Scale  Factor 

Hie  program  uses  the  following  definition  for  accelerometer 
scale  factors 


<v 

Where  a^  is  the  actual  specific  force  and  a^  is  the  accelerometer  output* 
i.o.*  an  increase  in  accelerometer  scale  factor  decreases  the  magnitude 
of  the  accelerometer  output.  The  scale  factor  consists  of  a linear  tern 
plus  a term  proportional  to  die  specific  force. 


$-19 


S =>  [1  + SO  + SlaJ 

“ [1  + SO  + SI  $0  )] 
P o 


The  scale-factor  coefficients  may  have  different  values  for  positive  and 
negative  values  of  specific  force.  The  scale- factor  equations  ares 


S 


S 


l*S0(+)-Sl(+)(K)  'oi° 

1 + s°(_)  * si(_,  (|  e0)  eo<0 


where 


S° <+>'  S°(-) 
S1(+)'  S1(-) 


positive  and  negative  linear  scale  factors, 
positive  and  negative  rate  sensitive  scale  factors. 


(9)  Velocity  Quantization 

The  pulse  generator  computes  the  number  of  output  pulses  produced 
by  the  output  angle  0q  plus  the  stored  angle  residual. 

0 “ 0 (n  + 1)  + 0_(n)  0_ (o)  « 0 

o R R 

If  (8  > 0)  S « 1 + S0(+)  + Sl(+) (|  0) 

If  (8  < 0)  S « 1 *1*  S0(-)  + Sl(-)  (|  0} 

n0  Inte9Qr  ^ 

? ■ qSn 

o 

0.,  (n  + 1)  » 0 • ? 

R 0 
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where 


% 

Oq  = quantized  indication  of  float  output  angle. 

= float  output  angle. 

0R  = quantized  angle  residual. 

S = scale  factor. 

q = A0q  quantization  level. 

n . = number  of  A0  pulses. 

0 o 

The  output  pulses  are  then  converted  into  an  indication  of  the 
change  in  velocity  (integral  of  specific  force)  along  the  input  axis  of 
the  accelerometer.  The  output  angle  is  a direct  indication  of  specific 
force. 


a. 

i 


0 


Thus,  the  integral  of  the  specific  force  along  the  accelerometer  input 
axis  over  the  simulation  time  step  is 

% 

<5v.  = a, At 

i i 

■ ?qn64t 

The  accumulated  incremental  velocity  is 


Av^  - E5v^ 

The  navigation  equations  reset  Av^  to  zero  each  time  the  incremental 
velocity  is  transformed  through  the  platform  attitude  matrix. 
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6.1. 3 Acceleromater  Compensation 


The  accelerometer  compensation  module  (ACOMP)  accepts  incremental 
velocity  from  the  simulated  SDF  pendulous  accelerometers  and  produces 
compensated  incremental  velocity  in  body  coordinates.  The  following 
errors  are  compensated: 

• Accelerometer  input-axis  misalignment 

• Scale-factor  error  (positive  and  negative). 

• Scale-factor  acceleration  sensitivity  (positive  and  negative. 

• Accelerometer  bias. 

• K ^ g-squared  coefficients. 

• Anisoinertia  torque. 

• Output-axis  acceleration  torque. 

• Lever-arm  effects. 

A flow  diagram  of  the  compensation  module  is  shown  in  Figure  6-7. 
The  compensation  equations  for  each  platform  axis  are  summarized  below. 

(1)  Scale-Factor  Compensation 

The  incremental  velocities  from  each  accelerometer  are  compensated 
for  linear  and  acceleration  sensitive  scale  factors.  The  accelerometer 
output  expressed  at  the  acceleration  level  is 

% -1 
a.  S (a,  + 6a) 

i i 

where  a^  is  the  actual  specific  force  and  ^a  is  the  additive  accelerom- 
eter error.  Thus,  the  estimated  specific  force  given  a^  and  estimates 
of  the  scale-factor  and  accelerometer  error  is 

ai  ” Sai  " 
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SCALE  FACTOR 
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Since  the  output  of  the  accelerometer  is  incremental  velocity  rather 
than  specific  force,  the  compensation  is 

= SAv^  - 6aAt 


Separate  scale  factors  for  positive  and  negative  specific  force  are 
provided  in  the  compensation  module. 


S 


S 


Av. 

1+S0<+>  +Sl(  + )(— }.  Av>0 


1 + SO 


Av. 

l. 


->  + sl<->(-ir>'  A < 0 


where 

SO^+j,  SO ^ j * positive  and  negative  linear  scale  factors. 

Sl^+j,  Sl^_j  » positive  and  negative  acceleration-sensitive 
scale  factors. 


(2) 


g-Sguared  Sensitive  Error  Compensation 


Compensation  is  provided  for  the  coefficients,  i.e.,  the 

terms  which  are  sensitive  to  the  square  of  acceleration  along  the  input 

axis.  The  g-sensitive  Kq  and  coefficients  are  not  compensated  nor 

are  the  K.  , K.  , K g-squared  terms.  The  acceleration  error  produced 
xo  xs  ss 

by  K±i  is 


<5a. 


The  g-squared  compensation  is 

6v  - K, .Av^/At 
g il  i 
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where  K is  the  compensation  coefficient  and  1/At  is  the  rate  at  which 
the  incremental  velocities  are  transformed  through  the  platform  attitude 
matrix. 

(3)  Acceleration  Bias  Compensation 

Acceleration  bias  produces  a constant  acceleration  error. 

6a.  = -B 

1 


The  bias  compensation  is 


6v.  = BAt 

l 

(4)  Anisoinei  tia  Torque  Compensation 

Anisoirertia  torque  produces  the  acceleration  error 


PUP 


I )W  U) 
x i p 


where 


I. 

i 


/ 


P 

I. 

i 


accelerometer  pendulosity. 

principal  inertias  of  the  float  about  (i,  p) . 


The  compensation  is 

5v.  = -^(1  - I . ) A0 . A0  /At 

i p p lip 

where  A0 . and  A0  are  incremental  rotations  about  the  accelerometer 
i P 

input  and  pendulous  axes. 
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(5)  Output-Axis  Acceleration  Torque  Compensation 

Output-axis  acceleration  produces  the  acceleration  error 

I 


where 

P » accelerometer  pendulosity. 

IQ  « principal  inertia  about  the  float  output  axis. 

The  compensation  is 

I 

6vi  = ~ [A0q ( n)  - 'A0o(n  - 1)]/At 
where  A0  is  the  incremental  rotation  about  the  output  axis. 


(6)  Lever-Arm  Compensation 

For  a rigid  lever  arm,  the  specific  force  at  the  location  of  the 
accelerometer  is  equal  to  the  specific  force  at  the  reference  position 
plus  the  centripal  and  tangential  acceleration  produced  by  vehicle 
angular  motion. 


b b 

' \ + ^ib  X 


. b b.  . . b b 
(u> . . x r ) +5)..  xr 

— ib  — —ib  — 


where 

a^  ■ specific  force  at  the  accelerometer  in  body  coordinates, 
b 

a^  ••  specific  force  at  the  eg  in  body  coordinates, 
r - lever  arm  in  body  coordinates. 
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The  objective  is  to  compute  the  specific  force  at  the  vehicle 

i.U 

center  of  gravity.  The  lever-arm  correction  for  the  i accelerometer 
is  obtained  by  evaluating  the  itb  component  of 


<$v£  = - [A0_g  x (A0_g  x £b)  + A20_b  x rb]/At 

where 

A26_g  - A£g(n)  - A6_g  (n  - 1) 

A8g  = incremental  rotations  from  the  simulated  gyros 

q 

The  gyro  outputs  A0_  are  the  values  before  gyro  compensation. 

(7)  Total  Accelerometer  Compensation 

The  compensated  incremental  velocity  along  the  accelerometer 
input  axis  is 


Av.  (comp)  = S(Av, ) + (ov  + Sv  + 6v  + 6v  + 5v  ) 
i i B g i oa  Li 


where 


<5v 


Sv 


<5v. 


<5v 


OA 


Sv, 


scale- factor  compensation, 
bias  compensation. 

g-squared  compensation. 

anisoinertia  torque  compensation. 
OA  acceleration. 

lever-arm  compensation. 


The  three  compensated  accelerometer  outputs  are  then  transformed 
from  the  accelerometer  input  axes  to  body  coordinates. 
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where 


. b _b.  a 

Av  a C iv 
— a — 


compensated  accelerometer  outputs  along  the  three 
accelerometer  input  axes. 

(nonorthogonal)  transformation  from  the  accelerometer 
input  axes  to  body  coordinates. 

The  trans formation  matrix  C*3  is  labeled  QMIS  in  the  program. 

a 

The  elements  of  this  matrix  must  be  consistent  with  the  rotation 
matrices  which  specify  the  (i,  o,  p)  axes  of  each  accelerometer  relative 
to  body  coordinates.  These  rotation  matrices  are  labeled  QABXr  QABY,  QABZ. 
All  of  the  matrices  are  obtained  from  the  initialization  data  set  for 
the  condensation  module. 
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SECTION  7 


VELOCITY-ATTITUDE  ALGORITHM  FOR  A 
5TRAPDOWN  INERTIAL  NAVIGATOR 


The  velocity-attitude  algorithm  is  the  second  or  middle  module  of 
the  three  main  modules  which  comprise  the  navigation  (and  attitude)  soft- 
ware for  a SD  INS.  The  inputs  to  the  velocity-attitude  (V-A)  module  are 
ideally  the  "incremental  velocities"  of  the  "center"  of  the  vehicle  in 
the  body  frame  over  a computation  cycle,  rAn)  , and  the  "incremental 
angles"  through  which  the  body  has  turned  with  respect  to  the  inertial 
frame,  expressed  in  the  body  frame,  over  a fast  cycle,  A0“b(n).  The 
inputs  from  the  inertial  sensors  are  the  raw  pulse  counts  which  are 
corrected  in  the  preceding  compensation  modules  to  give  the  computer's 

best  estimate  of  AV^n)  and  AO^.  (n) . 

ib 

The  functions  of  the  VA  module  in  the  INSS  are  threefold: 

a.  The  initial  value  of  the  transformation  from  the  body 
frame  to  the  inertial  frame,  C^(Q),  is  computed  using  the 
"loaded"  values  of  roll,  pitch,  and  yaw,  to  form  C^,  and 
the  "loaded"  values  of  wander  angle,  latitude,  and  longi- 
tude, to  form  c\  (The  "loaded  values  of  the  above  param- 
eters may  include  initial  errors  in  each.)  In  normal  sys- 
tem operation  C^(0)  would  be  computed  by  the  initial  align- 
ment routine. 

b.  The  strapdown  alignment  matrix,  C^(t),  is  updated  every 
computation  cycle  on  the  basis  of  the  incremental  angles, 
Ad“b,  from  the  gyros.  The  method  may  involve  either  a 
d.e.ra.  or  quaternion  update  (of  the  first,  second,  or 
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third  order)  but  in  any  case  the  updated  transformation 
is  expressed  as  a d.c.m.  for  use  in  the  AV  transformation 
and  for  attitude  computation. 

c.  The  incremental  velocity,  AV*5,  is  transformed  from  the  body 

frame  to  the  inertial  frame  every  computation  cycle,  by 

premultiplying  by  cj,  In  an  earlier  version  of  the  program, 

the  direction  cosine  matrix  per  se  was  not  generated  - since 

attitude  was  not  computed  - and  the  transformation  was 

accomplished  using  the  rotation  quaternion,  qj*,  and  its  con- 
i’* 15 

jugate,  qfe  , viz: 


AV* 


(7-1) 


One  additional  operation  is  performed  in  the  V-A  algorithm,  simply  for 
convenience  in  computing  attitude;  this  is  the  generation  of  the  trans- 
formation from  the  inertial  (i)  frame  to  the  "new"  earth  fixed  (e') 

9 * 

frame,  . This  transformation  accounts  for  the  multiply  defined 

inertial  frames,  i and  i*,  and  permits  the  direct  use  of  the  earth- 

c * 

fixed  to  "new"  LVWA  transformation,  C , , in  the  LVWA  navigation  algo- 

s 

rithm  to  complete  the  computation  of  the  "attitude"  matrix. 


7.1  Initial  Sfcrapdown  "Alignment"  Matrix,  cj^(O),  4,143  Initial 

. i 

Alignment  Quaternion,  q^tO) 

The  initial  body-to-iuertial  transformation  is  formed  as  the 

c i 

product  of  two  rotation  matrices,  and  Cc.  In  the  notation  of  Appen- 
dix B,  c£  is 


<S 


Z(it/2)X(ir)3<<J»  )Y(Q)X($)X(if) 
w 

" CP’SY  i -SR*SP*SY  - CR*CY  | -CR*SP*SY  ♦ SR*CY 
' t 

CP*CY  • -SR*SP»CY  + CP’SY  » -CR’SP’CY  - SR’SY 

» 

SP  1 SR* CP  ^ CR’CP 


(7-2) 
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where 


R = 

4> 

= roll  | 

f loaded 

F = 

e 

= pitch /values  may 

1 include  errors 

Y = 

w 

= yaw  J 

and  s = sin  and  c = cos 


In  the  V-A  module,  is  denoted  bv  QPB.  (A  lot  of  frames  are 

b 

the  P-frame  in  the  actual  INSS  program.) 

C1  at  t = 0,  and  hence  w.  t = 0,  is 
c xe 


C1(0)  = Y(fc)X(L)Z(a)  = Ce(0) 

c c 


CLO'CW 

- SLO«SLA*SW 

-CLOSW  - SLO»SLA*CW 

1 SLOCLA 

1 

1 SLA 

1 

=3 

CLA*SW 

CLA-CW 

-SLO*CW 

- CLO*SLA*SW 

SLO.SW  - CLQ»SLA»CW 

1 

| CLO’CLA 

(7-3) 


where 


LQ  “ i “ longitude 

LA  s L « latitude 

W “ a * wander  angle 


loaded 
values  may 
include  errors 


In  the  V-A  module,  C^(0)is  denoted  by  QIP.  Finally,  the  body-to- 

ic 

inertial  transformation,  formed  by  multiplying  the  above  two 

rotation  matrices,  i.o., 


C>> 


;{o)c^(o) 


(7-4) 


This  product  is  denoted  by  QIB  in  the  V-A  module 

In  the  case  where  a quaternion  update  algorithm  is  to  be  employed, 
it  is  necessary  to  convert  the  initial  body-to-inertial  d.c.m.  (or  rota- 
tion matrix),  C1  {0} , to  the  corresponding  unit  (or  rotation)  quaternion. 
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q^(0).  Later  (in  the  V-A  module),  when  a quaternion  update  is  performed, 
the  updated,  normalized  (unit)  quaternion  is  converted  to  a rotation 
matrix  for  AV-trans formation  and  attitude  computation.  The  link  between 
a unit  or  rotation  quaternion  and  the  corresponding  rotation  matrix  is 
the  interpretation  of  the  rotation  matrix  in  terms  of  the  direction 
cosines  of  the  axis  of  rotation,  and  the  angle  of  rotation  developed  in 
Appendix  C.  Calling  the  initial  alignment  matrix,  R = [rii] , the  initial 
(t  = 0)  quaternion  elements  are  computed  as  indicated  below; 


qQ(0)  = 

/ 11 

' ' 22 

s/ 

4 

<^(0)  = 

(^(Q)  - 

r23(0))/4q0(0) 

q2(o)  0 

l*13<°)  " 

r31(Q))/4q0(0) 

q3(0)  0 

(r21(0)  - 

rl2(0)1/4qO {0) 

’33' 


(7-5) 


7.2  Updating  the  SD  Alignment  Matrix  or  Quaternion 

Updating  C*(t)  or  q*(t)  from  t to  t + At,  requires  the  integration 
of  one  or  tho  other  of  the  following  equations; 


* i 1 b 

li  " Cb^ib  a update 


(7-6) 


or 


q*  * — for  a quaternion  update 


(7-7) 

where  the  superscript  dot  denotes  differentiation  with  respect  to  time. 


Obviously  the  simplest  fona  of  integration  which  ©ay  be  employed 
(first  order)  is  to  multiply  the  derivations  by  At  and  add  them  to  tho 
old  values,  viz; 
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7.2.1  First-Order  PCM  Update 


- 


c*(t  + at)  * c£(t>  + 


= c£[i  + a^b(t)at] 


(7-8) 


where 
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compensation 
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Then  the  C^(t  + At)  resulting  from  the  first-order  utxiate  is 


V 


t + At)  ~ 


fc  , c, , c,  . *1 

r i -A0n  *i 

11  12  13 

2 

C , CL..  c„. 

A0  1 -A0j 

21  22  23 
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(7-Ba) 


Unless  all  the  60^^  are  zero,  C^(t  + At)  will  not  be  a true  rotation 
matrix.  The  criterion  for  a rotation  matrix  is 
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I 


(7-9) 


or 


*v  —2  'v 

cc  » cc  * cc 


c'1c 


det  C = det  C * det  (CC)  = +1 


To  attain  (or  approach)  the  above  condition  the  updated  d.c.m.  is  (at 
least  periodically)  "orthormalized"  - to  the  first  order  - by  using  the 
following  algorithms 


^orfch 


1 'v 

C - - C (CC  - I) 


(7-10) 


7.2.2  First-Order  Quaternion  Update 

By  analogy  with  Eq.  (7-8),  we  may  write  the  first  order  quaternion 
update  as 

q£(t  + At)  S 


While  the  matrix  by  matrix  multiplication  in  Eq.  (7-8a)  is  well  defined, 
the  corresponding  quaternion  by  quaternion  multiplication  indicated  in 
Eq.  (7-11)  has  not  yet  been  defined. 

Let  us  generalize  the  expression  in  the  square  brackets  in  Eq.  (7-11) , 

jj 

[I  + At],  and  call  it  the  quaternion  p,  where* 

p “ pQ  + £ (7-12) 

A " * * 

and  £ a p^i  + Pjj  + 

AAA 

and  i,  j,  k are  the  frame  vectors  (unit  magnitude,  mutually  orthogonal, 

AAA 

right  handed)  (in  the  particular  case  i,  j,  k are  along  the  x,  y,  z 
axes  of  the  body  frame). 


q£(t)  + q£(t)  t 


ihti 


1 + (t)  At 

z 


(7-11) 
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A good  concise  summary  of  the  elements  of  quaternion  algebra  is 
presented  in  Section  2,2,1  of  Reference  3.  The  product  of  two  quater- 
nions, p and  q,  may  be  written  as 

= q0PQ  - SL*E  + <30E  + 3P0  + E^E  (7-1 3a) 

The  preceding  vector-scalar  form  may  be  handy  analytically  but  the  most 
efficient  form  for  mechanization  in  a digital  computer  is  as  the  product 
of  a 4 x 4 matrix  by  a 4 x 1 vector,  i.e.. 


qp  = 


The  purist  may  object  that  Eq.  (7-13b)  is  an  expression  for  the  elements 
of  the  quaternion  product  rather  than  the  product  per  se  which,  of  course, 
is  correct.  To  be  rigorous  then  the  left  hand  side  of  Eq.  (2-16)  should 
be  expressed  as 

j(qp)0  (qp)x  (qp)2  (qp)3 

The  rearranged  product  pq  is  given  below: 

Pq  = P0q0  - E • £L  + P0q  + £q0  + (7-14a) 

where  it  is  seen  that  the  sole  difference  is  in  the  sign  of  the  vector 
(cross)  product  term  since  £X£  = -£X£ 
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- — 
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MW 

(7-14b) 


Equation  (7-13)  is 
order  case  of  Eq.  (7-11), 


perfectly  general, 
we  have 


When  applied  to  the  first 
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Note  that  p is  not  a unit  quaternion  unless  to  is  zero;  hence  the 

ib 

quaternion  product,  pq,  will  need  to  be  normalized  by  dividing  each 
element  of  the  product  by  the  square  root  of  the  norm,  N(pq),  where 


N (pq) 


(pq>Q  + (pq)J  + (pq>2  + <pq>3 


(7-16) 


The  first  order  d.c.m.  and  quaternion  updates  and  (ortho) 
normalization  have  been  presented  by  way  of  an  introduction  and  to  point 
up  the  similarities  and  differences  in  the  simplest  cases.  The  first, 
second,  and  third  order  updates  for  each  are  presented  in  the  following 
two  subsections. 
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7.3 


Direction  Cosine  Matrix  Updating 


The  equations  for  higher  order  updating-based  on  the  availability 
of  incremental  angles  at  discrete  times  are  derived  in  (3)  and  the 
results  are  presented  here. 

The  matrix  to  be  updated  from  t to  t + At  is  denoted  by  C(t) , 
while  the  update  matrix,  of  order  i,  where  i = 1,  2,  3,  is  M^(t,  At), 
while  the  resultant  updated  matrix  is  (t  + At)  or 

C.,(t  + At)  = C(t)Mi(t,  At)  (7-17) 

The  A0.  of  the  previous  sections  are  shortened  to  0^,  while  the  A0^  from 
the  previous  pass  (At  earlier)  are  shown  as  0£.  Further, 

e2  ^ 02  + 02  + 02  (7-18) 


The  first  order  d.c.m.  update  (as  before)  is 


C (t  + At)  = C(t)M1(t,  At) 


where 


M (t,  At) 


The  second  order  d.c.m.  update  is 


C2(t  + At)  = C(t)M2(t,  At) 


(7-19) 


(7-20) 
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where 


_ 2 2 
■(e2 + 9 

0ie2 

6193  ' 

M2(t,  At)  = M^t,  At)  + | 

6162 

-(023  + 8 [) 

9283 

- 0193 

62e3 

-(8X2  ♦ 02)„ 

Note  that  the  second  order  update  approaches  R (see  Appendix  C)  as 

sin  6 , „ 1 - cos  0 . 1 

— - — ->  1 and r -*•  -r  . 


The  third  order  d.c.m.  update  is 


C3(t  + At)  = C(t)M3(t,  At) 
where  M3(t,  At)  « M2(t,  ^ 


(7-21) 


12 


I 2 * 

1 20V  + e.e, 

I 3 12 


- 6.6 


2 1 


-™\  ♦ 8l02*  - 0,6; 


20202  + e3e;  - el8;  I-<m\  * 0,0; 


(20262  * 03e;  - v; 


2026l  + 0,0;  - 03s; 


S3°2>  ! 


* * • 

In  the  third  order  update,  the  terms  (0^  - fl.,0^)  correspond  to  (w  x w) . 
It  was  found  (5)  empirically  that  the  inclusion  of  those  terms  in 
regimes  of  very  high  angular  acceleration ,u,  degraded  the  algorithm 
performance.  Accordingly  these  terms  may  be  omitted  at  the  user's  dis- 
cretion. (they  are  presently  "commented"  out  of  the  third  order  d.c.m. 
update.) 


The  program  flow  for  the  d.c.m.  update  utilizes  Eqs.  (7-19),  (7-20), 
and  (7-21)  in  order,  as  required,  followed  by  Eq.  (7-17),  and  ortho- 
normalization using  Eq.  (7-10). 
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7.4  Quaternion  Updating  and  PCM  Generation 

Regardless  of  the  order  of  the  update,  the  form  follows  Eq.  (7-13), 
while  the  equations  are  derived  in  (5).  The  quaternion  to  be  updated 
from  t to  t + At  is  denoted  by  q(t),  while  the  updating  quaternion,  of 
order  i,  where  i - 1,  2,  3,  is  p^  (t.  At)  and  the  updated  quaternion 
is  q^  (t.  At)  or 


q^  (t  + At)  = q(t)p^  (t,  At) 


(7-22) 


The  elements  of  p^  (t,  At)  are  all  functions  of  the  A0./2,  shortened  to 

^ * 

0./2,  while  the  A6./2  from  the  previous  pass  are  shown  as  0./2.  Simi- 

1 2 1 1 
larly,  0 from  the  d.c.m.  case  becomes 


(7-23) 


The  first  order  update  quaternion  (as  before)  is 


.(1) 


(t,  At)  = 1 


+ & 


(7-24) 


or 


(1) 


(1) 


(1) 


tl) 


The  second  order  update  quaternion  is 


p!2)(t,  At)  » pU)(t,  At) 


P(1)  (t,  At) 
7-11 


1 k . ^ 

2 \2  1) 

HD 


(7-25a) 


Since  the  second  order  portion  of  the  update  is  a scalar  it  affects 


only  pQ,  so 


(2) 

1 

0 

_ n(i)  i /e\2  _ i M2 

P0  2 \2/  ~ 1 2 \2/ 

(2) 

- n(1)  i - 1 ? , 

— p.  , i - 1,2,3 

(7-25b) 

The  third  order  update  quaternion  is 
p(3)(t,  At)  = p(2)(t.  At) 


(7-26a) 


Since  the  third  order 

update 

ss 

P(2) 
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>i31 
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( 2 ) 

P1  ’ 

(3) 

P2 
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p->- 

(3) 

Cm 

43)- 

is  a vector, 


it  does  not  affect  p , 

o 


so 


(7-26b) 


The  similarity  between  the  third  order  terms  of  the  update  for  the  d.c.m. 
and  quaternion  should  be  noted.  The  main  difference  is  the  use  of  0^/2 
in  the  quaternion,  whereas  0^  appears  in  the  d.c.m. 

Obviously,  the  simplest  way  to  formulate  the  algorithm  to  imple- 
ment Eq.  ( 7— 13b)  would  be  to  form  the  four  elements  of  the  update  quat- 
ernion, p^,  of  order  i - i.e.,  start  with  first  order,  add  second  and 
third  order  terms  as  required,  form  the  4x4  matrix  from  the  elements 
of  p^,  and  perform  the  indicated  matrix  by  vector  {(4x4)  x (4x1)] 
multiplication. 
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The  V-A  module,  uses  three  separate  inline-coded,  matrix  by  vector 
roul^iplication  routines  - one  for  each  order  of  quaternion  update.  How- 
ever, the  algorithm  used  is  correct,  even  though  it  takes  about  twice  as 
many  statements  as  the  simpler  and  more  elegant  algorithm  (see  reference 
(5) , in  subroutine  HSINTG) . 

After  forming  the  updated  quaternion,  p^(t.  At),  it  is  periodi- 
cally normalized  by  dividing  each  element  by  the  square  root  of  the  norm 
as  indicated  in  Eq.  (7-16)  that  is 


a = ^ 


q0  + ql  + q2  + q3 


I 


1/2 


^^nnrm  ~ ^.j/^*  i — 0,  1,  2,  3 

norm  i 


(7-27) 


where  the  q's  are  elements  of  the  updated  but  unnormalized  quaternion, 

and  q.  are  the  elements  of  the  updated,  normalized  quaternion, 

norm 

Next,  the  quaternion  is  converted  to  a d.c.m. , as  shown  in 

i (i) 

Appendix  c,  to  obtain  (t  + At).  Omitting  the  subscripts,  super- 

norm 

scripts,  argument  and  comments,  the  d.c.m.  is 


(7-28) 

This  completes  the  portions  of  the  V-A  module  which  depend  on  the 
particular  kind  of  update. 


*In  subsequent  sections,  it  is  more  convenient  to  use  "t"  as  the 
present  module  operating  time,  rather  than  (t  + At). 
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7.5  Incremental  Velocity  Transformation 

Once  C*  has  been  obtained  by  either  of  the  foregoing  means,  it  is 
employed  to  transform  the  incremental  velocities  - accumulated  from  (t  - At)  to 
t - from  the  body  frame  to  the  inertial  frame.  The  equation  which  is 
actually  mechanized  is 

Av^t)  = c£(t  - ^|)Avb(t)  (7-29) 

where  Av*5  (t)  is  given  by 

Avb(t)  4 f ab(i)dT  (7-30) 

“t  - At 

The  approximate  mid-computation  cycle  value  of  the  body-to-inertial 
transformation,  6^(t  - -j ) , is  obtained  simply  by  averaging  the  updated 
with  the  "old'1  value  from  the  previous  computation  cycle  i.e. 

* c£(t  - 4 i[cb(t)  + cb(t  ’ At)]  (7“31) 

While  the  approximation  in  (7-31)  is  only  first  order  (in  A9) , its 
effects  are  not  cumulative. 

The  final  operation  performed  on  the  Ay’s  before  output  is  to 
transform  them  from  the  i-frams  to  the  i'  frame  (the  inertial  frame  used 
by  the  navigation  algorithm),  and  to  sum  them,  if  the  navigation  algorithm 
is  operating  on  a longer  computation  cycle.  The  transformation  of  Av's  is 

Ay1*  <t)  * cJWtt)  (7-32) 


*This  transformation  is  denoted  by  DCMMID  in  program  mnemonics. 


where 


i * C*  *\j 
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(7-32b) 
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7.6  Body  to  "New"  Earth-Fixed  Frame  Transformation, 

For  attitude  computation  the  updated  must  be  transformed  to 

the  "new"  earth  fixed  or  e' -frame,  used  in  the  navigation  algorithm.  It 

6 ^ 

was  shown,  in  Appendix  B,  that  this  added  transformation,  Cj,  is  given  by 


Z(a)iet)S(Tt/2)X(ir/2) 

'swiet  0 

“ie*" 

c(i).  t 0 
ie 

“swiet 

0 1 

» 

0 

(7-33) 


where  the  value  of  t used  in  Eq.  (7-33)  corresponds  to  (t  + At)  in  the 
notaticn  of  this  section. 

o * 

The  transformation  which  is  output  to  the  navigation  module,  C , 

b 

or  DCM  in  the  V-A  mnemonics  is 


* 


See  Appendix  B. 
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SECTION  8 


8.  LOCAL  LEVEL  NAVIGATION  ALGORITHM  AND  ALTITUDE  COMPUTATION  FOR  A 
STRAPDOWN  INERTIAL  NAVIGATOR 


The  navigation/attitude  algorithm  is  the  third  or  final  module  of 
the  three  main  modules  which  comprise  the  navigation  (and  attitude) 
software  for  a SD  INS. 

The  function  of  this  module  is  to  compute  the  velocity,  position, 
and  attitude  (roll,  pitch  and  heading)  of  the  vehicle  at  the  end  of 
each  "slow"  computation  cycle  — on  the  order  of  2 to  20  Hz. 

To  accomplish  this,  the  module  is  initialized  at  the  start  of 
navigation  (t  = 0)  on  the  basis  of  velocity  and  position  data  supplied 
by  the  trajectory  module. 

During  each  succeeding  (non-initialization)  pass,  the  module 
accepts  incremental  velocities,  in  an  inertial  reference  frame,  summed 
over  the  slow  cycle,  from  the  velocity  and  attitude  algorithm,  and  the 
barometric  altitude  from  the  altimeter  module;  this  is  all  the  informa- 
tion required  for  velocity  and  position  computation. 

For  attitude  computation,  the  module  accepts  the  direction  cosine 
matrix  relating  the  body  frame  to  the  "new"  earth-fixed  frame  from  the 
velocity  and  attitude  module  and  combines  it  with  the  DCM  relating  the 
"new"  earth-fixed  frame  to  the  computational  frame  (the  "position" 
matrix)  and  extracts  the  computed  roll,  pitch  and  yaw.  The  latter  is 
corrected  by  the  wander  angle  (azimuth  of  the  computational  frame)  to 
obtain  the  heading  of  the  vehicle. 

8.1  Analytical  Development  of  the  Local  Vertical  Wander  Azimuth  Navi- 
gation and  Attitude  Algorithms 

The  actual  form  in  which  the  equations  for  the  computation  of  the 
earth-referenced  velocity,  position  and  attitude  of  a terrestrial  iner- 
tial navigation  are  mechanized  in  an  airborne  digital  computer  has  been 
the  object  of  considerable  study  for  the  past  quarter  of  a century. 
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A survey  of  the  algorithms  employed  in  several  aircraft  INS's  and 
the  algorithms  recommended  for  future  use  (in  tactical  aircraft)  was 
covered  in  sane  detail  in  (R-977)  . The  model  for  the  earth  and  its 
gravitational  field  is  that  of  the  WG2-72  ellipsoidal  earth  (R-977) ; many 
of  the  equations  relevant  to  this  model  are  presented  in  PROFGEN » which 
is  the  mission  profile  generator  that  ultimately  will  be  used  to  drive 
the  simulator. 

In  terms  of  the  first  reference  above,  the  navigation  and  attitude 
computations  used  in  the  INSS  most  nearly  approximate  those  that  would 
be  employed  with  a space  stable  inertial  measurement  unit  (SS  XMU) , 
i.e.,  AV's  are  accumulated  in  an  inertial  frame,  with  a local  vertical 
wander  azimuth  (LVWA)  computational  frame,  i.e.,  the  equations  of 
motion  are  actually  integrated  in  the  LVWA  frame.  Still  in  terms  of 
the  first  reference,  the  algorithms  used  in  the  INSS  correspond  to  the 
upgraded  algorithm.  It  should  be  noted  that  the  sign  of  the  wander 
angle  used  in  the  navigation  algorithm  has  been  reversed  with  respect 
to  the  reference  to  conform  with  the  remainder  of  the  program. 


With  these  preliminaries,  the  equations  mechanized  in  tho  naviga- 
tion and  attitude  computations  are  presented.  Note  that  the  first  sec- 
tion is  extracted  diroctly  from  (R-977) . 

8.2  Navigation  Equations  In  Local  Vertical  Frame 

Tho  fundamental  equation  of  inortial  navigation  for  a navigator 
employing  an  earth  relative  computational  framo  (c-spaco)  is  given  by 

vc  « f°  + gc  - (ftc  + 2 ft?  ) vc  (8-0) 

- - _ ec  is  - 

c 

where  v is  the  earth  relative  velocity  in  computational  frame 

* CJ  i c 

v is  the  time  rate  of  change  of  v 


ec 


ft? 

le 


is  the  specific  force  (the  portion  of  inertial  accelera- 
tion sensed  by  the  accelerometers)  expressed  in  computa- 
tional frame  coordinates 

c o c c c 

A G -ft,  ft.  r , is  the  resultant  of  mass  attraction,  G , 

- — ip  le  — — 

c c c 

and  centripetal  force,  -ft^g  £ » expressed  in  computa- 
tional frame  coordiantes 

is  the  earth-fixed  frame  which  is  rotating  at  a constant 
angular  rate,  u)^e,  with  respect  to  the  inertial  frame 

is  the  earth-centered  inertial  frame  which  is  regarded 
as  non  rotating  with  respect  to  the  "fixed"  stars 

is  the  skew-symmetric  form  of  the  angular  velocity  of  the 

computational  frame  with  respect  to  the  earth-fixed  frame, 

. c 

expressed  in  computational  frame  coordinates,  . 

is  the  skew- symmetric  form  of  the  angular  velocity  of  the 
earth-fixed  frame  with  respect  to  the  inertial  frame. 


expressed  in  computational  frame  coordiantes,  w 


—re 


The 


magnitude  of  w ? , w . , is  a constant  but  the  components 
3 — re  le 


vary  as  a function  of  latitude  (and  wander  angle) 


The  function  of  the  navigation  computer  may  be  considered  to  be  the 
integration  of  equation  (8-0) , to  obtain  the  earth  relative  velocity, 
vC,  followed  by  a second  integration  to  obtain  the  earth  relative  posi- 
tion. Since  the  earth  relative  position  is  expressed  mainly  in  angular 
coordinates,  latitude,  longitude,  (and  wander  angle),  it  is  necessary  to 
convert  the  level  components  of  linear  velocity  from  the  first  integra- 
tion into  angular  velocity  form  by  dividing  by  the  appropriate  variable , 

radii  of  curvature  to  obtain  w which  is  then  converted  to  the  skew- 

— ec 

symmetric  form,  ftC  . The  differential  equation  which  is  integrated  to 

GO 


obtain  the  angular  position  is  £3 


ce  ftc 

c ec 


(8-0a) 
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The  actual  latitude,  longitude  (and  wander  angle)  are  extracted 

from  C°  via  the  appropriate  inverse  trigonometric  functiono.  No  primes 
c 

(')  have  been  added  to  the  c,  e,  i frames,  since  usually  one  each  com- 
putational, earth-fixed,  and  inertial  frame  in  sufficient.  However,  the 
primes  are  employed  in  the  following  presentation  to  differentiate  the 
"new"  c,  e,  and  i frames  from  the  "old"  c,  o,  and  i frames. 

The  remaining  position  variable,  altitude  is  obtained  by  integrat- 

inq  the  vertical  component  of  v . An  external  altitude  reference  is 

required  to  eliminate  the  inheront  instability  of  any  pure  inertial 

vertical  channel.  Thus,  the  minimum  number  of  scalar  integrations 

required  to  mechanize  a local  vertical  wander  azimuth  navigator  is 

e 

thirteen  (13).  However,  only  five  of  the  nine  elements  of  Cc  are 
employed  directly  in  the  navigation  computations  so  frequently  only  six 
of  the  nine  elements  are  computed. 

8 . 3 Navigation  Initialization 

On  the  initialization  pass,  at  t =*  0,  the  navigation  module  obtains 
the  initial  position  and  velocity  from  the  trajectory  module , while  the 
initial  position  errors  are  supplied  from  the  module’s  own  initializa- 
tion file. 

The  initial  value  of  the  (LVWA)  computational  frame,  c',  to  "new" 

c * 

earth-fixed  frame,  o',  transformation,  C , , or  A,  is  computed  from 


e 

C , » A 

c' 

“ 2(H)  Y(L)  X(a) 

cLcK. 

! -cots.2.  - sasLcJl 

[ sasJl  - casLcfl. 

a 

cLs?, 

1 cotc&  - sasLstf, 

• -saeft  - casLsil 

sL 

1 

. sacL 

I cacL 

where  L is  the  initial  latitude  plus  errors 

i is  the  initial  longitude  plus  errors 

a is  tho  initial  wander  azimuth 


8-4 


This  matrix  is  the  repository  of  the  (horizontal)  position  infor- 
mation for  the  navigation  algorithm.  It  is  developed  in  Appendix  B. 


8.4  Initial  Velocity 

The  initial  earth -relative  velocities  in  an  ENU  coordinate  frame, 

£ 

v , are  confuted  with  their  initial  errors  and  transformed  to  the  LVWA 
computational  frame , c ' . 

I l l 

V = V + £V 

-nav  — — 


c'  c'  J l 

v - C0  v 
— nav  x-  -nav 


(8-2) 


where  C 


c’ 


c'  V 
**  c. , Cl 


and  C 


and  C 


c* 

i' 

V 


* x(a) 


* 2(tt/2)  X(u/2) 


(See  Appendix  B) . 


The  form  in  which  equation  8-2  is  mechanized  is 


vc  = X(a)  2(tt/2)  2(tt/2) 
-nav  L 


v 

— nav_ 


(8-2a) 


since  the  quantity  if  square  brackets  merely  respresents  a rearrange- 


ment of  the  elements  of  v 

—nav 


c' 

8.5  Initial  Angular  Rate, 

The  initial  radii  of  curvature  are  coirputed  using  the  exact  formu- 
lae and  the  parameters  of  the  WGS72  ellipsoid 


rp  - h + 
nav 


f , 2 2 \ 1/2 

C1_G  A3l) 


(8-3) 
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and 


(8-4) 


where 


r 

e 


is  the  earth's  equatorial  radius 

is  the  square  of  the  eccentricity  of  the  meridional 
ellipse 


A^=3L  is  the  sine  of  latitude 

h =h(0)+eh  is  the  loaded  value  of  the  altitude  above  sea  level 
nav 

r^  is  the  radius  of  curvature  of  the  prime  vertical  ellipse 

r is  the  radius  of  curvature  of  the  meridional  ellipse 

m 

The  angular  velocities  about  east  and  north,  to  , to  , are  computed 

£ N 


from 


and 


(8-5) 

(8-6) 


where  in  general. 


and 


c*  c* 

v « v-  cot  - v sa 

e i J 

c * c1 

v„  = v sa  + v,  ca 

N 2 i 


The  angular  velocities  about  the  horizontal  axes  of  the  con^uta- 

tional  frame  are  computed  by  transforming  the  w„  and  to  from  the  SL ' 

E N 

frame  to  the  c'  frame  f H 


0 


c A 
w , , = p 
— e c — 


X(a) 


to. 


to. 


N 


L toy,  A 0 


(8-7) 


where  the  p ^ are  the  components  of  the  angular  velocity  of  the  con$>uta- 
tional  frame  in  the  new  "earth-fixed  frame". 

l 

Equations  8-3  through  8-7  are  mechanized  in  a subroutine  (ANGVEL)- 
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in  the  compu- 


The  components  of  the  earth's  sidereal  rato, 
tational  frame  are 


9 


Ci'  Vie 


= Ceie30)ie 

nT 


e Cr 
3 c* 


0) 


ie 


e3  A 


w 


ie 


(8-8) 


T 

where  e^  is  a unit  vector  whose  elements  are  [0,  0,  1]  . 

c * 

All  equation  8-8  says  is  the  elements  of  are  formed  from 

the  elements  of  the  third  row  of  the  A-matrix  and  the  earth's  sidereal 

rate,  to.  , which  latter  is  directed  along  the  Z axis  of  the  i'  and  e'- 
ie 

frames . 


The  total  angular  velocity  of  the  "new"  computational  ('c)  frame 
with  respect  to  the  "new"  inertial  (i')  frame,  in  the  "new"  computa- 
tional frame,  0).°.’  , , is 
—l  c 


c * c ' 

0).  , , a W,  , , + 

—i  'c ' —i ' e ' 


c 

w. , . 

— i'c' 


(8-9) 
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and  are  the  gyro  torquing  rates  for  a gimballed  LVWA  IMU.  In  the 

present  SD  case,  these  become  the  "matrix”  torquing  signals  to  update 
c 1 

the  C , for  the  Av  trantormation . 

o * 


where 


Coriolis  and  centripetal  compensation  requires,  u ^ , 

C * Q* 

w = w , , + 2to . , , 

-cor  — e'c*  — x’e' 


(8-10) 


= w, , , + wr.  , 
— i ' c ' — i ’ e ' 


c 1 

and  the  form  of  the  correction  is  w X v or  (see  equation  (8-0) 

—cor  — 


C * c 9 

(u>xv)=  a)  Xv  = Q v 
— — -cor  — cor  — 


(8-11) 


0 i -oj 

| COr3 


W 


cor. 


-co  | 0) 

l_  cor2  | cor] 


0) 


cor„ 

i. 


-w 


cor. 


c' 

'l 

c' 

'2 

c' 


Since  the  correction  is  applied  at  the  velocity  level  all  the 
elements  are  multiplied  by  the  computation  cycle,  At. 

Equations  8-8  through  8-11  are  mechanized  in  a subroutine  (TORCOR) . 


\ 
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Several  other  computations  are  required  in  preparation  for  the  first 
normal  pass  through  the  navigation  algorithm. 


First,  the  computational  frame  to  earth  fixed  transformation, 

computed  at  t=0,  A(0) , extrapolated  to  the  middle  of  the  first  computational 

cycle,  using  the  initial  values  of  the  computational  frame.  A(0) ; or 
c * 

w , , (0) , and  a second  under  update. 

— e c 


= A(0)  ji+fl°!c.(0) 

2 

fr)  + 

- 

where 

0 -p3(0)  p2 (0) 

0 

P3(0)  0 0 

P2(0) 

-P2(0)  0 0 

_p3(0)_ 

A “..g.tW 


(8-12) 


and  At  is  the  navigation  computation  cycle. 

The  expression  in  the  curly  braces  is  computed  in  subroutine  AUP,  which 
accepts  as  inputs  0 and  0 where  (in  this  case) 

y z 

ey  = p2io)i| 

e2  = p3(o>% 

and  the  update  matrix 

I + 0 + i 02 

Ct 


r /0  2 + 0 

i - ( y.  z 

1 9y 

' 2 

0 

Le  2 

! e e 

z 

1 JL_ 

1 y z 

-0 

. 2 

10  0„ 

1 1-0  2 

U y 

| Y z 

| JL_ 

1 

i 2 

Based  on  the  initial  values  of  latitude  east  velocity  and  angle,  the  latter 
is  extrapolated  to  the  middle  of  the  first  computation  cycle  using  a first  order 
algorithm. 
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(8-14) 


sin 


cos 


[a(^] 

[“  m] 


- sin 


cos 


[a(0)]  + 6 Jjcx  cos  [a(0)] 

[a(0)J  - ~\  sin  [<*(0)] 


where  6 [a(%)  A a i~j)  ~ Ct(O)  = 3(0) 


At 

2 


Vfi(0)  sin  L (0)  At 

r cos  L(0) 

E 


Mid-computation  cycle  altitude  is  obtained  by  extrapolation  of  the  initial 
value  and  the  initial  vertical  velocity 


h(0)  vu(0)  ^ 


(8-15) 


Both  the  extrapolated  A-matrix,  A and  the  extrapolated  altitude, 
are  employed  in  the  computation  of  gravity  at  A~.  t in  t)  subroutine,  GRAV. 
Two  components  exist  - the  north-south  and  the  vertical  components.  The  east- 
west  component  is  zero.  However,  the  level  axes  of  the  computational  frame 
are  not  generally  north-south  and  east-west  and  must  be  resolved  through  the 
wander  angle  (or  azimuth  of  the  computational  frame) . 

The  vertical  (up)  component  of  gravity  is 

"[9o  + \ ,in2'£(^)  + % sl"4  £(t)] 

s(^)  * \ E2(%) 


(8-l6a) 


where 


sin2  £( 


The  level  component,  and  9 y which  are  east  and  north  for  a®0  are 
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$2  (t)  ' % fi  ^1) sin  £(t)cos  £ff)  si»  s(^f) 


(8- 16b) 


where 


and 


m- 

gn 

fi  (4) sin  £ 

A 

/At\  . A / 

At\ 

cos 

L 

(-)  sin  a ( 

2/  ~ ' 

A 

Mt\  A / 

At\ 

cos 

L 

(t)  cos  a l 

2/  = ‘ 

go 

= 

32.0877057 

ft/sec2 

= 

0.16939081 

ft/sec2 

-4 

2 

gL, 

7.5281x10 

ft/sec 

-8 

-1 

X 

s 

9.6227x10 

ft 

1 

-10 

-1 

gh. 

= 

6.4089x10 

ft 

2 

-15 

-2 

h3 

6.8512x10 

ft 

gn 

= 

1.63xl0~8 

ft2/sec 

extrapolated  A matrix  is  also  1 

(8-16c 


M4i) 


corrections  (call  TORCOR)  more  nearly  appropriate  to  the  mid-computation 
cycle. 

Finally,  the  incremental  velocities  in  the  "new"  inertial  frame, 
Av1(0)  are  set  to  zero,  i.e., 

Avi(0)  ■ 0 


These  are  DVI  in  program  mnemonics.  The  initialization  of  the 
navigation  algorithm  is  complete. 
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8.6  LVWA  Navigation  and  Attitude  Computations  for  Normal  Computation 

Cycle 

The  following  operations  are  performed  every  navigation  and  atti- 
tude) computation  cycle  following  initialization.  While  At  is  used  for 
this  module  as  well  as  all  the  others,  it  does  not  necessarily  follow 
that  all  the  At's  are  the  same.  What  is  necessary  is:  the  At  for 
"later"  modules  must  be  the  same  as  or  an  integral  multiple  of,  the 
computation  cycle  of  all  the  "earlier"  moudles  which  supply  the  input 
data  to  it. 

Note  that  in  the  navigation  algorithm  the  inertial,  earth-fixed, 
and  computational  frames  are  all  the  "new"  or  primed  frames,  unless 
stated  otherwise. 


1)  Transform  Av's  from  the  Inertial  to  the  Computational  Frame 
The  effect  of  the  two-step  process  is  to  mechanize 


AvC' (t)  = C?| (t-  Av1'  (t) 


or  Ave*  (t)  a cj,'(t-  Avi'(t) 


<8-17) 


where  C^, (t-  » Z 


At, 


and  Av  (t)  = CQ,(t-  ^ Av  <t) 


where  c°\  (t-  ~)  « A(t-  &} 

Gf  2 2 

^ At 

and  K(t-  — is  the  transpose  of  the  extrapolated  value  of  the  "new" 
computational  to  "new"  earth-fixed  frame  transformation  generated  during 
the  preceding  computation  cycle  (or  in  initialization).  Strictly 
speaking  it  should  be  denoted  as  £. 
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The.  components  of  gravitational  and  centripetal  force  at  (t-At) 
are  computed  using  the  WGS72  ellipsoidal  formula , and  the  values  of 
I c it  h. 

g? (t-At)  =-  [c  + C A3  (t-At)  + C A*  (t-At) 

1 1^0  g2  31  g4  31  J 

+ |l-[<^  - Ch  A^x (t-At) j h(t-At)  + 0^  h3 (t-At) 

92 (t-At)  » Ch  A31(t-At)  A32(t-At)  h(t-At) 


2)  Update  LVWA  Velocity 

The  Av,  Coriolis,  and  gravity  terms  are  used  to  updated  LVWA 
velocity  (per  equation  8-0) 

vc  (t)  = vc  (t-At)  + Avc  (t)  + w (t-  X vc  (t-At)  t + g°  (t- 
— — —cor  i — 

(8-18 


3)  Add  Vertical  Damping  Term  (to  Vertical  Velocity) 
c 


vf(t)  a v°'(t)  + [W-  ~)  + cvd2  6h(t-  —•)]  At 


-3  -2 

where  C , « 1,62x10  sec 

vdl 


Sa  (t-  —■)  » Cvd3  J*  6h(t)dT 


6h(t-  —•)  » hg(t-  ^|)  - fi(t-  ^|) 


At, 


At, 


(8-19 


Sa  and  Ah  are  extrapolated  values  from  the  previous  computation  cycle 


4)  Interpolated  and  Extrapolate  Computational  Frame  Velocity 

The  velocity  extrapolated  to  the  middle  of  the  next  comp,  cycle  is 

f’(t-  % = | vC*  (t)  - j vC\t-At)  (8-20) 

The  interpolated  velocity  is 

vC  (t-  ~)  = | [ vC  (t)  + vC  (t-At)j  (8-21) 

5)  Compute  Angular  Velocity,  f 

rThe  angular  rates  of  the  computational  frame  with  respect  to  the 

earth-fixed  frame  expressed  in  the  computational  frame  coordinates  at  the 

c * At 

mid  point  of  the  comp,  cycle,  to  , , (t -) , are  required  to  update  the  A 

e • - e c 2 

matrix,  C ,,  from  (t-At)  to  t. 
c 

The  following  inputs  are  supplied  to  the  angular  velocity  subroutine, 
ANGVEL: 


c'  At, 
v 


h (t- 


sin  a (t-  ^~) 


cos  a (t-  ^|) 


, 2r  , At, 
sin  L (t**  — j) 


All  except  the  first  are  the  extrapolated  values  from  the  previous 
cycle,  while  the  first  is  the  interpolated  value  from  the  present  cycle. 
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The  subroutine  computes  f (t j)  using  Equations  8-3  through  8-7. 

6)  Update  Transformation  from  Computational  to  Earth  Fixed  Frame 
(from  t-At)  to  t) . 

The  values  of  P(t-  —■)  are  converted  to  angular  form  by  multiplying 
by  At/  i.e: 

Px(t-  ^|)  At  - 0 
p2(t-  &>  At  = ey 
p3(t.  % At  - 9a 

The  second  order  update  matrix  is  formed  by  subroutine  AUP,  using  0y  and 
0^  above  as  arguments 

G * 

The  updated  A-matrix,  C , (t)  is 

c 

A(t)  = C®|(t)  - (t-At)  [i  + 9 f “ ej2  (8-22) 

where  the  expression  in  square  braces  is  defined  by  Equation  8-13. 

(2 ) 

7)  Orthonormalizo  A (t) 

(2) 

The  second  order  updated  DCM,  A (t)  / Is  orthonorwalised  every 
NORTH  computation  cycle. 

An(t)  » A{2)(t)  - ~A(2)(t)  [**>(t>  A^^t)  - ij  (8-25) 

The  orthonormalized  form*  A^(t)/  is  restored  as  A(t).  ■ 
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In  (6)  the  A-raatrix  was  updated  to  t,  A(t) , using  the  angular 
velocities  appropriate  to  (t-  ^r)  , P(t-  ~) . Now,  A(t)  is  extrapolated 

k . £ & 

using  the  same  P(t-  =j)  and  the  same  subroutine,  AUP,  with  arguments 


At 

2 


At.  At 
2 2 


0 

y 


P3(t- 


At  At 
2'  2 


0 

z 


The  extrapolated  A-matrix  is 


A(t  + -|j  - c®',  (t  + -~)  o c®!  (t)  fi  + o + \ ej2 


(0-24) 


where 


0 


f o -e.  9 . 1 


9)  Compute  altitude  and  Perform  Vertical  Damping  Calculations 


The  altitude  is  updated  not  only  to  account  for  the  mean  vertical 

A C * l\\~. 

velocity  over  the  computation  cycle,  (t-  ~) , but  also  to  include  a 

term  proportional  to  the  difference  between  the  external  and  inertially  - 

At 

derived  altitude  over  the  same  period,  <5h(t —) . The  equation  mechanized 

is 

h(t)  = h(t-A)  + £0°‘<t-  ~)  + Cvd  3h(t-  —)]  At  (8-25) 

The  altitude  is  extrapolated  to  the  midpoint  of  the  next  comp,  cycle  by 

ft(t  + —•)  « h(t)  + v°‘(t)  (8-26) 

Similarly,  the  external  (barometric)  altitude,  h , is  extrapolated  to  the 

O 

same  time  using 

V»  ♦ % " | »0(t)  - | hg(t-At)  (0-27) 

At  a A t 

The  predicted  altitude  error  at  (t  + —~l , 6h(t  + ~j)  is 

$h(t  + « fifl(t  ♦ ~)  - £(t  + (8-28) 

An  approximate,  scaled  (by  Cvd3)  integral  of  the  altitude  error,  So, 
is  computed  and  applied  at  the  acceleration  (actually  incremental  velocity) 
level  when  a third  order  loop  is  employed.  Its  effect  is  to  provide 
steady  state  compensation  for  vertical  accelerometer  or  altimeter  bias  errors. 
Xt  is  mechanized  as 

fa(t  + ^|)  - Sa(t-  Sh(t  ♦ —-)  At  (8-29) 

It  is  applied  to  the  vertical  velocity  as  indicated  in  6-i9. 
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10}  Extract  Latitude,  longitude  and  Wander  Angle  from  A Matrix 
The  procedure  is  outlined  in  Appendix  B. 


cos  L(t) 


L 


All(t)  * A2l(t) 


1/2 


(8-30) 


If  cos  L(t)  * 0,  latitude  is  +90°  and  longitude  and  wander  angle  are 
indeterminate  and  previous  values  are  output. 

f-  . , , w 


* L(t)  « tan 


-1 


Jl(t)  « tan 


-1 


A31(t) 

cosL(t) 


A21(t) 

Au(t) 


(0-31) 


(8-32) 


a(t)  » tan 


32  (t) 


A33(t) 


(8-33) 


During  the  succeeding  operations  and  the  next  computational  cycle,  the 
extrapolated  values  of  the  sine  and  cosine  of  the  wander  angle  will  be 
required.  They  are  obtained  from 


6a  A a(t)  - a(t-At) 
sa(t  + ~)  “ sa(t)  + ~ ca(t) 
cS(t  + ^~)  ■ coi(t)  + sa(t) 


(8-34) 


1 2 ) Transform  Computational  Frame  Velocities  to  END  Frame 
The  equation  to  be  mechanized  is 

v^ (t)  » c£, (t)  vc  (t)  (0-35) 

„ ^,<t)  cj;<t)  vc'(t) 

H c 

* 

Normal  navigation  outputs. 


This  is  the  inverse  of  equations  8-2  and  8-2a.  Both  transforraar- 
tions  are  given  in  Appendix  B,  hence 


VA(t) 

ta 

X (tt/2  ) 

2 (tt/2) 

X(a)  v° 

(t) 

r 

— 

m 
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_ 

V 

9 tO 
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ca 

-sa 

c 1 
V1 

(t) 
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(t) 

a 

0 

sa 

ca 

c1 

V2 

(t) 

V 

(t) 

l 

0 

0 

c • 

LV3 

(t) 

The  value  of  a above  is  a(t) . 


13)  Compute  Gravity,  Extrapolated  to  Mid-Computation  Cycle 


The  components  of  gravity  in  the  computational  frame  at  the  mid- 

t 


c 1 At 

point  of  the  next  computational  cycle,  g (t  + — ^-) , are  required  for 


Subroutine  GRAV  is  entered  with 


the  velocity  update  on  the  next  pass, 
arguments  &(t  + — |)  and  fi(t  + , Equations  9-16  a,  b,  and  c are 

mechanized  returning  g°  (t  + ~~)  and  sin2  £(t  + —) , 


14)  Compute  Extrapolated  Regular  velocity 


The  angular  velocity  of  the  computational  frame  with  respect  to  the 
At 

earth  fixed  frame  ,£(t  + — ~) , is  computed  for  use  in  the  upgraded  estimates 
of  the  Coriolis  corrections,  which  follow. 


The  arguments  supplied  to  subroutine  ANGVEL  are 


aC* 

v (t  + -g) 


ca(t  ♦ ^|) 


sa(t  + ~) 


s?£(t  + ^~) 


K<t  + 
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ANGVEL  implements  Equations  8-3  through  8-7 


15)  Compute  Extrapolated  Coriolis  Corrections 

The  Coriolis  corrections  for  use  in  the  next  velocity  update  are 
computed  in  subroutine  TORCOR. 

The  arguments  are: 


A(t  + 


£(t  + 


(t 


The  subroutine  mechanizes  Equations  8-8  through  8-11,  and  returns 
(w  x vc  ) At  at  (t  + —) , where  w cor  is  given  by  Equation  8-10. 

""•COT  a ■" 

16)  Reset  Incremental  Velocity  Input 

i ® 

The  Av  (t)  received  from  the  velocity-attitude  algorithm  is  reset 
to  zero  and  passed  back  to  the  V-A  module. 

This  anticipates  the  case,  which  often  occurs  in  practice,  wherein 
the  velocity-attitude  algorithm  is  exercised  several  times  for  each  pass 
through  the  navigation  algorithm. 


17)  Compute  "Attitude”  Matrix  and  Extract  Attitude  Angles 


Using  the  DCM,  or  C^,  from  the  velocity  attitude  algorithm, and  X 

C * Q * 

or  Cgl ,from  the  navigation  algorithm, the  "attitude"  matrix,  is 
formed: 


c£(t)  = c°,(t)  c®(t) 


(8-37) 


Equation  8-37  is  developed  and  interpreted  in  Appendix  B,  where 
the  final  result  is 


Cfa  = vOr/2)  ZODw)  ?(6)  x (<f>)  X(tt) 


s0 


C0s4> 


s^i  c0  ( -c(J»  c<j)-siji  s0s<j> 

w f 

I cij;  c0  | c<j>-c^  s9s(j) 

W * W W 


C4  C5 


C0C({) 


C^wS(j)-S^wS0C(J) 

-st^s^-c^sOc^ 


(8-38) 


o DTEM  in  program  mnemonics 


From  the  above  (the  elements  of  the  first  row  and  column  of  ) 
the  "attitude"  angles  are  extracted,  viz. 


cosO 


(c 


i * 


(8-39) 


If  cosO  is  0 then  0 = +90°  and  i|)  and  $ are  indeterminate  and  the 

“ w 


previous  values  are  output.  Otherwise, 

K " tan'1  (J~) 

* $ « tan 


0 « tan 


~i(  ci ) 
\cos0  ) 


(8-40) 

(8-41) 

(8-42) 


Normal  navigation  outputs 
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Finally,  the  true  heading,  \p,  is  computed  from 

* = ij>w  - a (8-43) 

where  a is  obtained  from  equation  8-33  and  the  argument,  t,  has  been 

omitted  for  brevity  in  equations  8-38  through  8-43. 

* 

The  Av  (t)  received  from  the  velocity  attitude  algorithm  is  reset 
to  zero  and  passed  back  to  the  V-a  module.  The  module  operating  time 
is  stored  for  output,  then  incremented  for  the  next  pass. 

The  nine  stored  (*)  quantities  and  time  are  output  to  the  evalua- 
tion module,  and  the  navigation  module  awaits  the  next  pass. 


* 

Normal  navigation  outputs. 
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SECTION  9 


EVALUATION  OF  ERRORS 


The  evaluation  Module  (EVA)  is  the  final,  minor  module  of  the 
INSS.  Its  functions  are 

(a)  To  store  the  position,  velocity,  and  attitude  from  the 
trajectory  module  at  the  time  of  each  navigation  and 
attitude  computation  cycles. 

(b)  To  store  the  position,  velocity,  and  attitude  from  the 
navigation  module  after  every  navigation  computation  cycle. 

(c)  To  output  the  trajectory  data  (50  points  at  a time)  to  the 
trajectory  data  and  output  the  resulting  differences  or  errors 
(50  points  at  a time)  to  the  line  printer  and  plot  file. 

The  net  result  is  a time  history  of  the  nominal  trajectory,  i.e,, 
without  any  perturbations  from  the  environment  module,  and  an  instantaneous 
“error"  time  history.  As  presently  structured,  “errors"  include  the 
effects  of  the  linear  and  angular  vibration  on  the  simulated  strapdown 
navigator,  but  not  on  tho  nominal  trajectory.  This  is  an  obvious  defi- 
ciency  of  the  present  program  since  the  total  effect  of  its  vibration  en- 
vironment should  bo  reflected  in  the  position,  velocity,  and  attitude, 
comprising  the  nominal  trajectory,  as  reported  by  EVL. 
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APPENDIX  A 


ROTATION  MATRICES 


A. I Introduction 

A simple,  compact,  unambiguous  notation  describing  transformation 
from  one  orthogonal  axis  frame  to  another  in  terms  of  one  or  more  single- 
axis rotations  is  presented  in  this  section. 

The  notation  used  in  describing  the  vectors  and  the  coordinates 

( 4 ) 

transformations  is  based  on  notation  used  by  Britting.  The  defini- 

tions and  the  notations  presented  in  this  Appendix  are  employed  in 
Appendix  B to  describe  the  coordinate  frames  and  transformations  occur- 
ring in  the  INSS. 


A. 2 Single-Axis  Rotations 

A coordinate  frame,  i,  is  defined  by  specifying  three  orthogonal 

A A A 

unit  vectors  X^,  Y^,  Z^,  which  form  a right-handed  system,  i.e., 

A A A A A AA  A A 

X,  H ¥.  a 2,,  Y,  X Z,  a X,  , z.  X X.  a Y, 

i i i i i i i i i 


X, 


« 0,  Ya 


Zi  a 0,  2i  • Xi  a 0 


where  "x"  denotes  the  vector  product  and  "•**  the  scalar  product  of  two 
vectors . 

The  **i"  is  the  shorthand  name  of  the  frame,  and  is  one  or  more 
numbers  or  letters,  or  some  combination  thereof.  The  practice  of  refer- 
ring to  the  unit  vectors  defining  a frame  as  i,  j,  k»  u,  v,  w;  M,  N,  Oj 


A-l 


AAA 

X,  Y,  Zj  etc.  will  not  be  used  because  there  are  too  many  frames  and  one 
would  soon  run  out  of  triples  of  letters  for  which  the  sequence  is  obvious. 
Further,  the  existence  of  an  X-axis,  a Y-axis,  and  a Z-axis  in  each  frame 
is  indispensable  to  the  concept  of  single-axis  rotations. 

Consider  two  Cartesian  coordinate  frames  with  a common  origin. 

Denote  the  first  frame  by  i and  the  second  by  j.  Then  the  unit  vectors 
defining  the  two  frames  are  X^,  Y^,  Zy  and  X^,  Y^ , Z^.  If  one  axis  of 
one  frame  coincides  with  the  corresponding  axis  of  the  other  frame,  i.e., 
if 

A A A A A A 

X,  ? X,  or  Y.  » Y.  or  Z.  = Z, 

i j i 3 i j 

then  the  transformation  from  one  frame  to  the  other  is  at  most  a single- 
axis rotation.  The  axis  which  is  common  to  the  two  frames  is  the  axis 
of  rotation  and  tne  angle  between  the  corresponding  non  cosanon  axes  is 
the  angle  of  rotation  or  argument.  The  sign  of  the  angle  of  rotation  is 
positive  if  the  non  common  axes  of  the  first  ("from")  frame  move  in  a 
clockwise  direction  to  bring  themselves  into  coincidence  with  the  cor- 
responding axes  of  the  second  ("to"  frame,  looking  in  the  positive  direc- 
tion along  the  common  axis  (right-hand  screw  rule) . 


(1)  X-axis  Rotation,  X (a) 

A A 

If  X.  ° X.r  the  transformation  is  an  X-axis  rotation*  and  if  the 

* aJ  a A a 

angle  between  Yi  and  Y^,  and  and  Zy  is  a and  if  a is  positive,  then 
transformation  from  i to  j is  written  as  X(a) . This  is  illustrated  in 
Figure  A-l. 

The  equations  relating  the  unit  vectors  of  the  j frame  axes  to  the 
unit  vectors  of  the  i frame  axes  are* 

A A 

Xj  - Xi 

AAA 

Y,.  ® cos  o Y^  + sin  a Z^ 

A A A 

Z ^ “ -sin  a Y^  + cos  a Z^ 
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or  in  matrix  form 


X. 
3 

Y. 
3 


L 3 J 


10  0 
0 cos  a sin  a 

0 -sin  a cos  a 


LZi 


x(a)  4 


1 

0 

0 


0 0 

cos  ot  sin  a 

-sin a cos  a 


A A 

Xj-*i 


Figure  a-1.  Rotation  about  X-axis. 


4 

In  this  particular  case,  the  roation  from  i to  j,  denoted  by  is  X(a). 
It  is  often  helpful  to  represent  a number  of  single-axis  rotations  in  a 
signal  flow  graph  form  as  shown  in  Figure  A-2. 


© s *©  - ©^  "s  © 

“from"  axis  of  “to" 
frame  rotation  frame 


Figure  A-2.  Rotation  about  X-axis  flow  diagram. 
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The  arrowhead  indicates  the  "to"  frame,  changing  the  direction  of 
the  arrow  corresponds  to  changing  the  sign  of  the  argument  or  transposing 
the  rotation  matrix. 


(2)  Y-axis  Rotation;  Y(8) 

Similarly,  if  = Y^  the  transformation  from  the  j frame  to  the  k 
frame  is  a Y-axis  rotation.  For  a positive  argument,  3 , this  is  written 
as  Y(0)  and  is  illustrated  in  Figure  A- 3. 


Y (0) 


A 

S3 


cos  3 


-sin  0 


Lsin  0 


cos  0 J 


Figure  A-3.  Rotation  about  Y-axis. 


Figure  A-4.  Rotation  about  Y-axis  flow  diagram. 


(3)  Z-axis  Rotation,  Z(  ) 

Finally,  if  Z.  = Z0,  the  transformation  from  the  k frame  to  the 
£ K * 

i,  frame,  c,  , for  a positive  argument,  y,  is  written  as  Z(y)  and  is 

A 

illustrated  in  Figure  A-5. 

cos  y sin  y 

Z(y)  ^ -sin  y cos  y 
. 0 0 


Figure  A-5.  Rotation  about  Z-axis 


Figure  A-6.  Rotation  about  Z-axis  flow  diagram. 


A.  3 The  General  Rotation 

If  the  three  single-axis  rotations  are  cascaded,  the  transformation 

l 

from  the  i frame  to  the  i frame,  C^,  may  be  written  as 

k j 

ci  * cj  ci 


Zly)  VIS)  X(a) 


Note,  that  the  subscripts  cancel  the  superscripts  from  lower  left 
to  upper  right,  leaving  only  the  lower  right  subscript  (“from")  and  the 
upper  left  superscript  ("to").  This  latter  notation  appears  in  Brittin^ 
The  vector  is  designated  by  a small  letter  with  one  or  more  subscripts, 
if  required.  The  frame  in  which  a vector  is  defined  is  denoted  by  a 
superscript  and  transforming  a vector  v from  the  i frame  to  the  & frame 
is  expressed  as 

ft  JL  i 
v » v 

For  a vector  with  notation  employing  two  subscripts,  such  as  the  angular 
velocity  of  the  x frame  with  respect  to  the  j * frame,  expressed  in  the 

4,t. 

k ‘ frame  is  expressed  by  uk  . . 

Several  of  these  general  rotation  matrices  occur  In  the  simulation 
and  mechanization  of  a strapdown  ins.  The  required  transformations  arc 
presented  in  Appendix  B. 


j ^vs.-v<  ^ Vs.  - 


A prerequisite  to  establishing  a simple,  consistent  set  of  single- 
axis rotations  relating  the  frames  of  interest  is  to  define  the  coordinate 
axes  of  one  frame,  and  then  relate  the  others  to  it  via  the  appropriate 
rotations . 


APPENDIX  B 


COORDINATE  FRAMES  AND  TRANSFORMATIONS 
AND  THEIR  USE  IN  THE  INSS 


B.  1 Introduction 

In  any  inertial  navigation  which  provides  position  and  velocity 
with  respect  to  the  earth  in  a local  geodetic  vertical  north  pointing 
(LVN)  frame  and  the  attitude  (roll,  pitch,  and  heading)  of  the  vehicle 
with  respect  to  the  same  LVN,  frame,  there  are  generally  four  fundamental 
coordinate  frames  of  interest,  namely « 

1.  The  Platform  Frame 

• Defined  by  the  input  axes  of  the  inertial  sensors. 

2.  The  Body  (or  Vehicle)  Frame 

• Defined  by  the  longitudal,  lateral  and  normal  axes  of 
the  vehicle. 

3.  The  Local  Vertical  North  Frame 

• Defined  by  the  local  (geodetic)  vertical,  east  (or  west) 
and  north  (or  south)  axes. 

• Fixed  with  respect  to  the  earth  at  a point  below  (or  above) 
the  vehicle. 
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4.  The  Inertial  Reference  Frame 


• Defined  by  the  earth's  polar  axis  and  the  equatorial 
plane  (at  the  start  of  navigation) . 

• Nonrotating  with  respect  to  the  celestial  sphere. 

Only  the  first  of  these  fundamental  coordinate  frames,  (1) , is  a 
function  of  the  configuration  and  mechanization  of  the  particular 
terrestrial  inertial  navigator.  In  fact,  the  three  classical  mechani- 
zations of  an  INS  are  based  on  the  platform  frame  coinciding  with  one 
of  the  three  remaining  fundamental  coordinate  frames.  Thus,  if  the  plat- 
form frame  coincides  with  the  body  frame,  we  have  a strapdown  system? 
if  the  platform  frame  coincides  with  the  local  vertical  north  frame, 
we  have  a local  vertical  system;  if  the  platform  frame  coincides  with 
an  inertial  frame,  we  have  a space  stable  system.  The  latter  two  types 
imply  stabilization  of  the  platform  or,  in  these  cases,  the  stable  element 
by  gimbals  and  control  loops  driven  by  gyro  and  other  signals;  the  strap- 
down  system  mounts  the  instrument  cluster  directly  to  the  body.  It  is 
this  type  of  inertial  measuring  unit  (IMU)  which  is  considered  in  this 
report. 


B.2  Digression 

An  added  complication  has  crept  into  the  mechanization  o£  strap- 

down  systems  due  to  the  fact  that  they  were  preceded  chronologically  by 

local  vertical  systems.  While  local  vertical  systems  work  very  well  over 

most  of  the  earth's  surface,  both  computational  and  physical  problems 

arise  in  the  vicinity  of  the  earth's  poles.  The  computational  problem 

occurs  in  the  calculation  of  the  longitude  rato  which  is  of  the  form 

V / ( ft  cos  L),  whore  V„  is  the  east  velocity  with  respect  to  the  earth, 

E 

R is  the  radius  of  the  earth,  and  cos  L is  the  cosine  of  the  latitude. 

As  the  latitude  approaches  90  degrees  (at  the  poles) , cos  L approaches 
zoro  and  the  longitude  computation  blows  up.  The  physical  limitation 


B-2 


arises  in  the  computation  of  the  azimuth  gyro  torquing  command,  i.e., 
the  rate  at  which  the  platform  must  be  rotated  about  a vertical  axis  to 
maintain  one  of  its  level  axes  north.  The  torquing  rate  is  of  the  form 
[(VE/R)  tan  L + ohQ  sin  L] , where  u^e  is  the  sidereal  earth’s  rate, 
and  tan  L and  sin  L are  the  tangent  and  sine  of  the  latitude.  Again, 
as  the  latitude  approaches  90  degrees,  the  torquing  rate  may  approach 
infinity.  To  circumvent  these  singularities,  the  local  vertical  wander 
azimuth  (LVWA)  mechanization  was  developed,  wherein  the  "azimuth"  gyro 
torquing  was  reduced  or  eliminated  by  incorporating  its  effects  in  a 
"wander"  angle  which  becomes  one  of  the  Euler  angles  in  a direction 
cosine  matrix  (d.c.m.).  Since  a d.c.m.  never  becomes  singular,  the  ex- 
traction of  latitude,  longitude,  and  wander  angle  from  the  Euler  angles 
of  the  d.c.m,  avoids  the  problem  of  the  north-slaved,  local  vertical 
mechanization.  However,  both  the  longitude  and  wander  angles  are  in- 
determinate at  the  poles. 

The  computational  and  physical  limitations  of  a north-slaved, 
local  vertical  system  need  not  arise  in  strapdown  (or  space  stable) 
system,  since  the  gyros  are  not  torqued  to  maintain  a preferred  orienta- 
tion with  respect  to  the  earth.  The  IMU  of  a strapdown  system  outputs 
incremental  velocities  and  angles — the  integrals  of  specific  force  and 
angular  velocity1 — in  the  body  (platform)  frame.  If  the  incremental 
velocity  outputs  are  transformed  to  the  inertial  frame  and  integrated 
using  Cartesian  coordinates,  the  position  and  velocity  may  be  expressed 
in  geographic  coordinates  without  resorting  to  the  introduction  of  a 
wander  angle. 

However,  in  the  present  simulation,  the  use  of  a LVWA  computational 
frame  was  specified  by  customer,  so  the  incremental  velocities  in  the 
inertial  frame  are  transformed  to  LVWA  computational  frame,  where  the 
navigation  algorithm  appears  the  same  as  it  would  for  a gimbaled,  local 
vertical  XHU- — except  that  the  gyro  torquing  signals  become  the  matrix 
torquing  signals  for  the  inortial-to-LVWA  transformation. 
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The  purpose  of  this  digression  has  been  the  explanation  (if  not 
justification)  of  the  introduction  of  an  additional  coordinate  frame, 
the  LVWA  computational  frame.  This  frame  coincides  with  the  LVN  frame 
when  the  wander  angle,  a,  is  zero.  When  the  nominally  "north"  axis 
(for  a = 0)  is  rotated  to  the  west  of  true  north,  the  wander  angle  is 
positive. 

B.3  The  Simplest  Case 

If  one  were  seeking  the  simplest  mechanization  of  a strapdown 
system  for  terrestrial  navigation,  there  would  be  only  three  coordinate 
frames  of  interest.  They  are  as  listed: 

1.  The  platform  (pj , end  body  (b)  (or  vehicle)  frame. 

2.  The  local  vertical  north  frame. 

3.  The  inertial  reference  frame. 

In  Figure  B-l,  a gross  "circle  diagram"  illustrates  the  three 
coordinate  frames  of  interest  (denoted  by  the  small  circles) , and  the 
intervening  transformation  (denoted  by  the  rectanyular  boxes) . The 
circle  diagram  is  a pictorial  representation  of  the  fact  that  the  product 
of  the  rotations  around  a closed  path  in  the  same  sense  is  identity  or 

b i • i* 

<V  <v  =D  ■ 1 (B_1) 

With  the  direction  of  the  rotations,  as  shown  by  the  arrowheads  in 
Figure  B-l,  which  corresponds  to  premultiplying  both  sides  of  Eq,  B-l 
by  C^, , the  attitudo  matrix,  we  have 

„1*  i'  JL*  _ 

eb'ci'°i'cb  • cb  1 IB'2> 
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B-l . Transformations  required  for  strapdown  (SD  INS 
compu -..Lions  - simplest  case). 


cc 


(B-3) 


.(4  ) 


is  used  throughout.  Any  orthogonal  trans- 


whore  Drifting 's  notation 
formation,  , corresponds  to  ono  or  tho  product  of  two  or  more  single- 
axis  rotations,  as  indicated  in  Appendix  A. 


To  maintain  the  simplost  relationships  among  the  fundamental 
coordinate  frames,  requires  that,  when  the  angles  specified  by  the  rota- 
tions are  all  zero,  the  coordinate  frames  are  coincident.  This  is 
accomplished  by  specifying  the  coordinate  axes  of  one  frame  and  relating 
it  to  the  others  by  tho  simplest  sequence  of  single-axis  rotations. 


If  the  Z axis  of  the  inertial  reference  frame  is  defined  to  lie 
along  the  earth's  polar  axis  (EPA)  with  the  positive  end  at  the  north 
pole,  and  the  X axis  is  defined  by  the  intersection  of  the  equatorial 
plane  and  the  meridian  plane  (observer's  or  Greenwich  at  time  t = 0) , 
then  the  Y axis  lies  in  the  equatorial  plane,  90  degrees  east  of  observer's 
or  Greenwich  at  t = 0.  This  defines  an  earth-centered  inertial  (ECI) 
reference  frame  for  all  t >_  0 , If  the  first  definition  is  used,  then  the 
longitude  computed  is  the  change  in  observer's  longitude  since  t = 0, 
and  the  observer's  initial  longitude  relative  to  the  Greenwich  meridian, 

ZQ,  is  simply  an  additive  constant.  Positive  longitude,  or  change  in 
longitude  , is  to  the  east  along  the  equator,  i.e.,  it  is  a positive 
rotation  about  the  EPA.  The  rotation  of  the  earth  with  respect  to  the 
ECI  frame,  is  about  the  same  axis  and  in  the  same  sense  as  longi- 

tude. 

Then,  to  define  the  location  of  a point  on  the  earth,  P,  the 
observer's  location  at  time  t,  a rotation  about  an  easterly  axis  parallel 
to  the  equatorial  plane,  through  the  geodetic  latitude,  L,  is  required. 
Since  latitude  is  defined  to  be  positive  in  the  northern  hemisphere, 
the  sense  of  the  rotation  is  negative  in  the  northern  hemisphere. 

In  the  compact  notation  of  Appendix  A,  the  transformation  from  the 

inertial  reference  frame  to  the  LVN  frame  (observer's  position  with 

1 ' 

respect  to  the  earth  at  time  t) , , , is  expressed  as  the  product  of 

two  single-axis  rotations: 

C*;  » Y(-L)ZU  + wiQt)  (B-4) 

The  X,  Y,  Z axes  of  the  LVN  (£')  frame  previously  defined  are  Up,  East, 
North,  respectively. 
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If  the  aircraft  were  headed  north  and  the  wings  and  nose  were 
level,  all  the  attitude  angles  would  be  zero  and  the  body  frame  would 
coincide  with  the  LVN  frame.  However,  if  the  aircraft  performs  a 
climbing,  right-hand  turn,  all  three  "attitude"  angles  (roll,  <{>? 
pitch,  6;  heading,  0)  will  assume  positive  values.  The  simplest 
transformation  from  the  UEN  (LVN)  frame  to  a body-fixed  frame  involves 
three  successive  single-axis  rotations  through  the  heading,  pitch,  and 
roll.  They  are  as  follows: 

1.  The  first  is  a negative  rotation  about  the  X or  "Up"  axis  of 
the  LVN  frame  through  the  heading  angle,  0;  i.e.,  when  0 goes 
positive,  the  nose  of  the  aircraft  is  rotating  from  north  to 
the  east.  This  rotation  is  denoted  by  X(-0)  = X(0)  , where 
the  tilde  (~)  denotes  the  transpose. 

2.  The  second  is  a positive  rotation  about  the  Y axis  of  the 
intermediate  frame  (which  is  displaced  from  the  UEN  frame 
by  x(-0) ; i.e.,  the  Y axis  is  East  if  0 » 0)  through  the 
pitch  angle,  0.  When  the  pitch  angle  goes  positive,  the 
nose  is  elevated  above  the  horizontal  plane.  This  rotation 
is  denoted  by  Y(9). 

3.  The  third  is  a positivo  rotation  about  the  forward-pointing, 
longitudinal  (Z)  axis  of  the  aircraft  through  the  roll 
angle,  0.  Roll  is  positive  when  the  right  wing  dips  below 
the  horizontal  piano.  This  rotation  is  denoted  by  2(0), 

b' 

The  complete  UEN (LVN)  to  body  transformation,  C^, , is  given  by 

cj|  » Z(0)Y (0)  X(-0)  (B-b) 

Note  that  the  X,  Y,  2,  "body"  axos  in  this  case  are  "Up"  through  the 
the  canopy,  and  through  the  right  wing,  and  forward  through  the  nose, 
respectively. 
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The  attitude  angles  as  previously  defined  in  Eg.  (B-5)  correspond 
to  those  obtained  from  perfect  angle  encoders  on  an  ideal,  three-axis, 
local  vertical,  north-slaved  IMU,  where  the  gimballing  order,  from  the 
outside  in  is  roll,  pitch,  and  azimuth.  This  gimbal  arrangement  is  used 
on  cruise  vehicles,  where  the  pitch  angle  does  not  approach  90  degrees. 

On  tactical  aircraft  with  all-attitude  capabilities,  a redundant,  inner- 
roll  gimbal  is  added  frequently  with  limited  freedom — which  drives  the 
outer-roll  gimbal  so  as  to  maintain  the  inner-roll  gimbal  orthogonal  to  the 
pitch  gimbal  (except  during  certain  maneuvers) . This  gimbal  arrangement 
and  the  resulting  "attitude"  matrix  is  illustrated  in  Figure  A-9  in 
Section  A. 2 of  Appendix  A,  Volume  XI. 

In  a strapdown  system,  the  gyros  are  usually  mounted  with  their 

input  axes  along  (or  parallel  to)  the  aircraft  body  axes.  The  gyros 

ideally  sense  the  angular  velocity  of  the  body  frame  with  respect  to 

the  inertial  frame,  expressed  in  the  body  frame.  The  outputs  of  the 

gyros  represent  the  integrals  of  these  angular  velocities  over  a short 

interval  of  time.  So  if  the  initial  orientation  of  the  body  with  respect 

i<  b' 

to  the  inertial  frame,  C^COfc  or  (^,(0^  inortial-to-body,  is  known 
and  one  properly  keeps  track  of  the  changes  in  orientation  with  time, 
then  C^, (t)  is  always  known. 

Now,  premultiplying  Eq.  (B-4)  by  Eq.  (B-5)  wo  have 


2<$)Y<e)X(-ii>)Y(~L)ZU  luiot> 


(B-6) 


1* 

Actually,  fcho  transpose  of  the  matrix  given  in  Eq.  (B-6)  or  C^,  is  com- 
X>utod  in  the  "velocity-attitude"  algorithm  which  computes  (or  can  compute) 
the  inertial  to  local  vertical  transformation,  C^,.  Then  by  premultiply- 
ing the  former  by  the  latter,  we  obtain  the  "attitude"  matrix,  i.c., 


Since  the  form  of  Eq.  (B-7)  is  known  from  Eq.  (B-5) , the  attitude  angles, 
i|»,  6,  # may  be  extracted,  as  shown  in  Section  A. 4. 4 of  Appendix  A, 

Volume  II.  This  simple  case  is  illustrated  in  Figure  B-2) . 

B.4  The  Actual  Case 

The  coordinate  frames  and  transformations  actually  employed  in 
the  INSS  are  considerably  more  complicated  than  the  simple  case  just 
presented.  The  complication  arises  due  to  the  introduction  of  seven 
additional  single-axis  rotations  required  to  form  the  body- to- computational 
frame  transformation,  . Two  of  these  additional  rotations  occur  due 
to  the  use  of  a LVWA  computational  frame,  with  the  introduction  of  the 
wander  angle,  a,  which  appears  twice  in  the  sequence,  while  the  remaining 
£tve  are  fixed  rotations  which  occur  in  ones  or  twos  to  accommodate  the 
various  definitions  of  local  level  or  inertial  frames  which  appear  in  the 
different  modules.  The  overall  nominal  body-to-computational  frame  circle 
diagram  is  shown  in  Figure  B-3. 

Before  discussing,  this  circle  diagram  in  detail  it  may  be  instruc- 
tive  to  examine  when  all  the  variable  anglos,  <j>,  0,  iji^,  a,  L,  i., 
and  “ie*  ate  aQ3CO*  rQmambering  that  in  the  simple  case  all  the  coordinate 
frames  would  be  coincident.  This  situation  is  depicted  in  Figure  B-4  where 
all  the  remaining  six  frames  are  regarded  as  LVN  framos.  Under  those 
circumstances  five  different  LVN  frames  appear  since  the  north,  west, 
up  framo  occurs  twice. 

c * 

In  Use  Figure  B-4,  , with  all  variable  angles  zero  is  given 

by 


all  variable 
angles  zero 


0 0.1 


0-10 


10  0 


Z(n) V(w/2> 
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C = z = e = i 


ENU 
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Figure  B-4.  C,  in  INSS  with  all  variable  angles  zero, 
b 

Hence,  it  requires  a total  of  seven  fixed,  single-axis  rotations  to 
complete  the  trip  around  the  circle  diagram  when  all  the  variable 
angles  are  zero.  With  a little  planning,  all  these  rotations  could 
have  been  eliminated,  as  illustrated  by  the  "Simplest  Case" — or  as 
used  in  Reference  2. 

In  practice,  it  is  common  to  encounter  different  definitions  of 
the  same  kind  of  coordinate  frames,  although  seldom  as  extreme  as 
shown  in  Figure  B-3. 
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B. 5 The  LVWA  (Computational) -to-Body  Frame  Transformation, 
c£  or  QBC 

This  transformation  relates  a wander  azimuth,  local  vertical, 
computational  frame  (corresponding  to  the  stable  element  of  a perfect 
LVWA  navigator)  to  the  body  frame.  It  is  used  in  the  interim  trajectory 
module  as  part  of  the  transformations  required  to  generate  the  incremental 
velocities  and  angles  which  would  be  output  by  an  ideal  strapdown  IMU. 

It  is  expressed  as  the  product  of  six  single-axis  rotations — three  fixed, 
and  three  variable: 

Cb  = X{  tt)  X{  A)  Y ( 6)  2 ( i{j  ) X ( ir)  Z ( rr/2) 
c T w 


where 

<f>  = the  roll  (Euler)  angle. 

0 = the  pitch  (Euler)  angle. 

iD  = the  yaw  (Euler)  angle  of  the  wander  azimuth  platform, 
w 

Performing  the  indicated  operations , we  have 
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A 

c = cos 
A 

s = sin 
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It  should  be  noted  that  3x3  matrices  are  handled  in  the  program  as  nine 
element  vector  (9x1) , where  the  e ants  are  stated  rowwise,  i.e. 


or 


Cc  or  QBC  = 


"°11 

C12 

C13_ 

’C1 

C2 

C3" 

C21 

C22 

C23 
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C4 

C5 

C6 

rH 

m 

u 

1 

C32 

C33 

_C7 

C8 

C9 

QBC  = {C1C2C3C4C5C6C7C8Cg} 


B . 6 The  Local  Vertical  North  (LVN)  to  Local  Vertical  Wander 
Azimuth  (LVWA)  Transformation,  or  QCP 

This  transformation  relates  an  ENU  (LVN)  frame  to  an  ideal  LVWA 
frame.  It  is  a single  axis  rotation  about  the  Up-axis  through  the 
wander  angle,  a,  which  is  positive  when  the  Y-axis  of  the  LVWA  frame 
is  west  of  north. 


C*  = Z (a)  = 


ca  sa  0 

-sa  ca  0 

0 0 1 


In  the  Interim  Trajectory  Module  (ITM)  the  incremental  velocities 
(integrals  of  specific  force)  are  generated  in  an  East,  North,  Up  (ENU) , 
Local  Vertical,  North-pointing  (LVN)  coordinate  frame,  then  are  trans- 
formed to  the  body  frame  (which  would  be  NWU  if  all  the  variable  angles 
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were  zero).  In  the  notation  of  the  circle  diagram,  Figure  B-3,  this 
transformation  is  denoted  by  C.  . However,  as  indicated  in  the  heading, 
the  program  mnemonic  is  QBP  (P  as  is  in  perfect  LVN) . 


b b c 

CZ  = Cc  °Z 


The  corresponding  body  frame  specific  force,  a , is  given  by 


-b 

a 


b -Z 
C Z 3 


B. 7 The  Earth  Fixed  to  LVN  Transformation,  C^-  or  QPE 

This  transformation  relates  an  earth-centered  earth-fixed  frame 
with  its  Z and  X axes  in  the  equatorial  plane,  and  its  Y-axis  along 
the  north  earth's  polar  axis  (EPA)  to  the  perfect  LVN  frame  with  its 
X,  Y,  and  Z axes  al< 
rotations,  namely, 
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B‘8  "The  Inertial”  to  Earth  Fixed  Transformation,  or  QEI 

This  is  the  single-axis  transformation  that  relates  "the  inertial" 
frame  to  the  earth-centered,  earth-fixed  frame.  It  is  a position  rota- 
tion about  the  north  EPA  through  an  angle  which  is  the  product  of  the 
earth's  sidereal  rate,  <^e,  and  the  time  in  navigation,  t.  It  is  given 
by 


Cl  = Z (ui.  t)  = 
1 ie 


CO),  t 

ie 


-SO),  t 
le 
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The  angular  rate  sensed  by  the  gyros  of  a strapdown  system  is 

that  of  the  body  with  respect  to  the  inertial  frame  expressed  in  the 
— b 

body  frame,  The  output  of  the  gyres  of  the  SD  IMU  is  the  integral 

of  the  aforementioned  angular  race,  over  some  time  interval.  At,  or 


(t)dt 


To  evaluate  this  integral,  the  product  of  all  the  four  rotation  matrices 
is  required;  i.e.,  the  transformation  from  the  inertial  frame  to  the 


b (tnl 

body  frame  at  t and  t , , C.  n 
1 n n-1  l 


C. 
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where  the  tilde  (~)  denotes  the  transpose,  then  the  product  of  the 
previous  matrices 


b(tn)  _ b(tn)  i 

b(t  i b(t  n 

n-1)  n-1) 

permits  the  evaluation  of  A0  , as  indicated  in  Appendix  C. 

n 

JL' 

D . 9 The  Inertial  to  "New"  Earth-Fixed  Transformation,  C or  Q 

In  the  earlier  discussion  of  the  various  LVN  frames,  which  could 
be  considered  to  be  implied  by  the  numerous  fixed  rotations,  it  was  not 
specifically  noted  that  some  of  these  frames  arose  due  to  two  different 
definitions  of  the  "inertial"  frame  and  the  "earth-fixed"  frame. 

The  velocity-attitude  algorithm  (module)  in  the  INSS  is  responsible 
for  updating  on  the  basis  of  the  A0's  received  from  the  gyro  compensa- 
tion module,  so  the  AV' s received  from  the  accelerometer  compensation 
module  may  be  transformed  from  the  body  frame  to  the  inertial  frame  used 
by  the  navigation  algorithm  (module) . In  addition,  the  velocity-attitude 
algorithm  must  pass  to  the  navigation  algorithm  the  body-to-"new"-earth- 

Q * 

fixed  transformation,  C,  , for  use  by  the  latter  in  the  attitude  compu- 

b 

tat ions . 

e ' . e ' 

The  additional  transformation  involved  xn  generating  is  , 
where  e'  is  the  "new"  earth-fixed  frame,  and  i'  is  the  "new"  inertial 
frame.  These  frames,  i and  e,  and  i’  and  e',  coincide  (in  pairs)  when 
t = 0.  The  first  pair  is  based  on  the  X,  Y,  Z axes  of  the  LVN  frame 
being  along  East,  North,  Up,  respectively,  while  the  second  case  assumes 
the  X,  Y,  Z axes  of  the  LVN  frame  are  along  Up,  East,  North,  respectively. 
The  two  cases  are  illustrated  in  Figure  B-5. 
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Figure  B-5. 
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The  transformation  from  the  "old"  inertial  frame,  i,  to  the  "new"  inertial 

i ' 

frame,  i',  is  C.  where 

C1'  = Z ( — tt/ 2 ) X ( — tt/2 ) 

1 

= Z (tt/2)  X(ir/2) 

From  the  "new"  inertial  frame  to  the  "new"  earth-fixed  frame  the  trans- 
0 1 # 

formation,  , , is  simply 

G * A 

c , = Z(W.  t) 

i ' ie 
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Finally,  the  required  transformation  from  the  inertial  to  "new"  earth- 


0 * 

fixed  frame,  , is  given  by 


1 0 * *1  * 

c.  = c , c! 

1 1 i 


Z (to . t)  Z ( tt/2 ) X (tt/2) 
le 


CO),  t 

SOJ . t 

0 

' 0 

-1 

cf 

" 1 

0 

0* 

le 

xe 

-SO),  t 

CO) . t 

0 

1 

0 

0 

0 

0 

-1 

lr 

le 

0 

0 

1_ 

_ 0 

0 

1_ 

_ 0 

1 

0_ 

SO),  t 

le 


CO),  t 

xe 


L o 


0 

1 


ca>.  t 
le 


-SO) . t. 

xe 


0 J 


B. 10  The  ("New")  Earth- Fixed  to  ("New")  Computational  (LVWA) 

Frame  Transformation,  Cc*  or  A^  or  A 
e ' 

This  transformation  (or  its  transpose)  is  loosely  and  generally 
referred  to  as  "the  direction  cosine  matrix",  or  d.c.m. , in  inertial 
navigation  systems  employing  a LVWA  navigation  algorithm  for  the  simple 
reason  that  it  may  be  the  only  rotation  matrix  or  d.c.m,  occurring  in 
the  navigation  software  for  a gimballed,  LVWA  INS. 

c 1 

The  elements  of  C^,  are  the  repository  for  the  earth-relative 

position  (latitude  and  longitude  and  wander  angle)  which  are  the  three 

Euler  angles  relating  the  "new"  earth-fixed  frame,  e',  to  the  "new" 

c * 

computational  frame,  C'.  C^,  is  given  by 
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X(a)Y(-L)ZU) 


c ~ A 
C . or  A = 
e' 


where 


= X(a)Y<L)ZU) 

10  0 cL  0 sL  Ci  si  0 
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a =>  the  wander  angle  positive  to  the  west  of  north 
L “ the  geodetic  latitude  positive  in  the  northern  hemisphere 
i =*  the  geodetic  longitude  positive  to  the  east  of  Greenwich 


From  the  elements  of  A,  position  and  wander  angle  are  obtained  from  the 
following  relationships 


a = tan 


-1  A32 


i,  = tan 


-1  A21 
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L = sin  1 (&31) 


2 2 2 2 

Note  that  a and  £ are  indeterminate  if  A, , + A„,  = A..  + A,,  = 0, 

11  21  32  33 

i.e.,  for  CL  = 0 or  L = +90°. 


B. 11  The  "Attitude"  Matrix  or  Body  to  "New"  Computational 

Frame  Transformation,  CJ-'or  DTEM 


This  transformation  relates  the  body  frame  to  the  "new"  LVWA 
computational  frame  employed  in  the  navigation  algorithm.  By  the 
appropriate  inverse  trigonometric  relationships,  the  "attitude"  of  the 
vehicle  (body)  is  extracted.  "Attitude"  included  roll,  pitch,  and  yaw. 

The  latter  must  be  corrected  by  the  wander  angle  to  provide  the  "true" 
heading. 

0 I 0 1 

In  the  INSS,  Cfa  is  formed  as  the  product  of  from  the  velocity- 
c * 

attitude  algorithm,  and  C , from  the  navigation  algorithm,  or 

0 


<£ 


I 


» 


This  looks  simple  enough  until  one  reexamines  the  circle  diagram, 

c * 

Figure  B-3,  and  notes  that  involves  16  single-axis  rotations  and 
8 variable  rotations  through  a,  L,  £,  t (each  of  the  latter  4 argu- 
ments appearing  twice) . One  would  expect  that  the  eight  aforementioned 
variable  rotations  would  "cancel"  out  in  pairs,  leaving  the  resultant 
purely  a function  of  the  three  attitude  angles,  6,  and  0,  with  suit- 
able rearrangement  due  to  the  fixed  rotations. 
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In  the  "simplest"  case,  shown  in  Figure  B-2,  one  can  see  by 
inspection  that  this  is  the  case,  but  in  Figure  B-3,  it  is  not  so  ob- 
vious that 


<£'  = F(\|»  , 0,  $) 

D W 

C * 

Letting  A = Z + m^t  and  tilde  (~)  denote  transpose  may  be  written 

as 


C 


c 

b 


where 

<f  » Z (ir/2)  X ( it ) Z ( xp  ) Y (8)  X(<j>)  X ( tt ) 
b w 

and 


c * • 

It  may  be  shown  by  expanding  Cc  , as  indicated  by  the  successive  brackets, 
that  only  the  fixed  rotations  survive,  and 

Cc'  = Z (ir/2)  X(ir/2) 
c 


Hence 


Cf'  = Z(TT/2)X(ir/2)Z(Tr/2)X(Tr)Z(^  ) Y (9)  X(<{>)  X(ir) 

b w 
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Expanding  the  four  fixed  rotations  which  are  adjacent  to  one  another,  C 
we  have 


Cp  = Z(TT/2)X(Tr/2)Z(TT/2)X(Tr)  = Y(tt/2) 

o 


c 1 

The  final  form  of  the  "attitude"  matrix,  C.  , is 
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o 


Roll,  pitch,  and  yaw  are  obtained  from  the  elements  of  the  first  row 
c * 

and  column  of  C£  , namely 

*w  - t“'1  (4) 


2 2 2 2 

Note  that  and  are  indeterminate  if  CL  + C,  =»  C.  + C,  » Oj  i.e.,  for 

w 2 3 4 7 

C0  ® 0 or  0 *=  ±90°. 

B.12  True  Heading 

It  may  be  shown  that  the  true  heading,  <!>,  is  the  difference  be- 
tween the  yaw,  and  the  wander  angle,  a,  i.e., 

to  « !|)  - a 

w 

In  the  mechanisation,  the  position  and  wander  angle  are  extracted  before 
the  "attitude"  angles,  so  the  current  value  of  a is  available  when  required. 
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D • 1 3 The  "Wander  Angle  Matrix",  or  "New"  LVWA  Computational 

Frame  to  LVN  (UEN)  Frame  Transformation.  C* 

c • 

This  single-axis  transformation  is  used  to  rotate  the  computational 
frame  velocity  to  an  LVN  frame.  It  is  given  by 

*10  0 * 

V 

C , = X(-a)  = X(a)  = 0 ca  -s a 

0 sa  ca 


and  its  uses  is 


l * 
v 


c 


V 


I 


B. 14  Other  Transformations 

Another  fixed-angle  transformation  is  employed  to  express  the  LVN 

velocity  in  the  ENU  frame,  when  it  is  given  in  the  UEN  frame.  The 

l 

transformation  is  , and  is  given  by 


cj,  = X (tt/2)  Z (tt/2) 


0 10 
0 0 1 
10  0 


and  its  use  is 
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This  latter  form  is,  of  course,  the  acme  of  simplicity  to  code, 
but  it  is  a frequent  source  of  trouble  in  interpretation. 

If  nothing  else,  the  presence  of  so  many  redundant  transforma- 
tions points  up  a very  real  problem  in  systems  integration,  where 
multiply-defined,  similar,  coordinate  frames  are  almost  certain  to 
occur.  The  construction  and  use  of  the  appropriate  circle  diagram(s) 
is  recommended  in  these  instances. 
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APPENDIX  C 


ROTATION  MATRIX  IN  TERMS  OF  AXIS  AND  ANGLE  OF  ROTATION 


A direction  cosine  matrix  (or  rotation  matrix)  is  unique,  its 
interpretation  is  variable.  As  the  first  name  implies,  the  direction 
cosine  matrix  is  a 3x3  array,  the  rows  or  columns  of  which  are  the 
direction  cosines  of  the  axes  of  one  orthogonal  frame  in  another . The 
rotation  matrix  may  arise  in  the  case  where  a vector  is  rotated  about 
a line,  i.e.,  a unit  vector  with  direction  cosines  in  the  "body”  frame 
of  n^,  n2,  n3,  through  an  angle  6. 

In  this  case,  the  rotation  matrix  from  the  new  body  frame  to  the 
old  (before  the  rotation)  body  frame,  R(6) , is  given  by  (3) 

2 

R = I + N sin  0 + n (1  - cos  0) 

where 


C-l 


Another  equivalent  expression  for  R is 


0 

~n3 

n2 

“ 2 

ni 

nin2 

nin3 

R s=  I COS  0 + 

n3 

0 

"ni 

sin  6 + 

nin2 

2 

n2 

n2n3 

”n2 

ni 

0 

‘nln3 

n2n3 

2 

n3 

— 

— 

_ 

•PM 

The  direction  cosines  of  the  axis  of  rotation  n^,  n2,  n3#  may  equally 
well  be  expressed  in  terns  of  the  components  of  the  angle  of  rotation, 
0,  along  the  axes  of  the  coordinate  frame,  as 


n,  ** 


whore 


/i  2 2 

/01  0*  + 0~ 


In  terms  os  6 and  its  components,  IS  and  tT  become 


02 


and 


and  R becomes 


2 2 
-(02  + 03> 


0ie2 


ei02 


2 2 

-(03  + e,) 


6193 


0203 


(02  + e3> 

1 ~ : U - cos  8) 


3 , °1°2 

— sin  0 * — c-  (1  - cos  8)  ■ 

• ' a • 


12 


0103 


0203 


-(0?  + 6?) 


0, 


sin  8 ♦ (1  - cos  8)  — sin  0 + (1  - cos  8) 


<S  j * *J>  ' 

1 (1  - COS  0) 


2 3 

— T"  i*  - t 


®2  ®3®1  - 81  M,  - (0  * Of) 

- 5-  sin  0 ♦ — 5-  (1  - po 8 0)  ~ sin  0 ♦ -£-i  (1  - cos  0)  1 - i.  (1  - cos  8) 

_ 9 8 02 
One  should  bear  in  mind  that  any  rotation  matrix  may  be  considered  to 

have  the  form  of  R (it  is  not  limited  to  small  valuos  of  8)  although 

the  simplifications  and  approximations  which  are  reasonable)  and  as  0 

2 

becomes  smaller,  are  evident.  If  terms  on  the  order  of -01  are 
negligible  then 


R 


02  + 0 


1 


since  1 


cos  0 + 


sin  0 


-*>  1 


9 


2 


6 


+ 1 


6 


In  general,  given  R,  the  angles  0^,  6^,  0^,  are  extracted  as  follows 


v - r 
32  31 


2 

r - 

r 

13 

31 

2 

r , - 

r. 

21 

12 

5 S^-n  Q 
31  0 


6 sin  0 
2 0 


a sin  0 
63'-¥~ 


Squaring  the  above  three  equations  and  adding  gives 


,„2  , a2  . a24  sin  0 . 2 A 

<0,  + 6„  + 0,)  -x — = sin  0 

1 2 3 0 2 


From  this  define 


sin  0 = + t4in2  0 


and  hence 


0 = sin  ^ (sin  0) 


If  an  inverse  sine  routine  is  not  available  or  convenient,  the  trace 
of  on  the  sum  of  the  diagonal  elements  is  used  to  define  cos  0, 
namely, 


r.  + rA<5  + r » TrR  = 1 + 2 cos  0 

XI  J J 


hence. 


cos  0 « (TrR  - l)/2 
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and  the  inverse  tangent  routine  may  be  used  to  evaluate  9.  Note  that 
there  is  no  loss  of  generality  in  assuming  the  positive  sign  for 

/ . 2 r 

/sm  9 , since  a change  of  sign  here  would  be  reflected  in  the  signs 
°f  G^,  02,  and  0^  effectively  flipping  the  axis  of  rotation  and  the 
sense  of  the  angle  of  rotation  simultaneously.  j 

This  interpretation  of  a rotation  matrix  is  used  to  generate  the 
incremental  angles  which  would  be  output  by  ideal  single-degree-of- 
freedom  gyros  in  a strapdown  system,  at  t where  the  values  of  C* 

are  qiven  at  t and  t - At.  In  reference  (3),  the  procedure  described 

n n b(n) 

above  is  used,  but  the  INSS  simply  forms  the  product  matrix  as 


b(n) 

b(o) 


cb(n) 


cb(o) 


T 


cb(n) 


~b{o) 


R 


where 


n -+  nAt,  0 -+•  (n  - l)At 

Then,  three  off-diagonal  elements  are  extracted  and  divided  by  At,  and 

th 

the  resultants  are  called  the  average  angular  rates  over  the  n com- 

2 

putation  cycle.  In  addition  to  the  scaling  error  of  6 /6,  this  method 

also  has  additive  errors  in  0.  of  9.9/2. 

l j K ' 

C.  1 Relationship  Between  a Rotation  Matrix  and  a Unit 
Quaternion 

The  general  rotation  matrix,  R,  shown  in  expanded  form  in 
Eq.  (3-  ),  may  equally  well  be  expressed  exlicitly  in  terms  of  the 
direction  cosines  of  the  axis  of  rotation,  n^,  n^,  and  n^,  and  of 
trigonometric  functions  of  the  half  angle  of  rotation,  9/2. 

Using  the  half-angle  trigonometric  identities 


3in  0 


- , 0 G 

2 sin  — cos  — 


=2 


2 s 


0. 

2 


1 


COS  0 a 


2 sin 


2 £ 
2 


2 s 


2 0_ 
2 


05 


1 


er  6 “ n!»  02/8  “ n2'  03/e  “ n3'  nl  + n2  n3 


and  substituting  the  above  into  the  expression  for  R,  we  have 


1 - (n2  + n2) 

«-as» 


e e. 

-n3(2  s j c J) 
+ nin2)2  s2  f) 


„ 6 0. 

n2(2  bjc-) 

z_  2 e, 

+ n3n  (2  s j) 


n3 12  s f 0 f 1 

+ “ln2(2  f1 


1 - (n2  + n2)  (2  s2  |-) 


-ni12  s f 0 l> 

+ n2n3(2  s2  f) 


-n2<2  s | o |) 
2 6, 

+ n3ftl(2  s 2) 


nx(2  sf  of) 

+ n2n3(2  *2  I1 


1 - (n2  + n2) (2  s 


Noting  from  Reference  (3)  that  the  four  elements  of  a unit  quaternion 
q#  may  be  expressed  as 


e e 

cos  y c J 


. e e 

nlsin2  ' "l82 


. e e 

"2  sln  I V n2  S 2 


, 9 0 

a n,  sin  r B n,  s •? 


|;j  i 

I!  | 

K I 


*l  I 

Id  I 

S'  ! it 

*■  ] It 

ijl  |! 


.A  I 

|j  | 


R may  be  written  in  terms  of  the  elements  of  the  unit  quaternion  as 


1 - 2(q2  + q2) 


2(qiq2  " q3q0}  2 (q3qi  + q2qQ) 


m i 

6!  * 


R"  2(q3q0  + qiq2) 


1 - 2(q2  + q2) 


2(q2q3  ‘ Vl* 


" V 


|y  1 1 
rf\  i 

II  !: 

W.  f 

I ; i 

k i 


§4.  I 

h l 

l l 

t I 


2%h  ' *2V 


2 (q2q3  + q0qi} 


1 - 2(q2  + q2) 


2 2 2 2 

Since  qQ  = 1 - q^  - q^  - q^,  ^or  a un^fc  quaternion,  the  elements  of 
the  quaternion  may  be  found  from  the  elements  of  the  corresponding  rota- 
tion matrix,  as  follows: 

1 + ru  * r22  * r33  ■ 1 + = 4(1  - *1  - 4 - «J>  * *4 


hence 


q^  ® /(I  + TrR)/4 


% “ (r32  ‘ r231/4  «0 


92  ’ (r13  " r31)/4  % 


<r21  " r12l/4  % 


Both  the  above  conversions  are  employed  in  the  Velocity  Attitude  algorithm 
when  a quaternion  update  algorithm  is  used. 
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APPENDIX  D 


INERTIAL  COMPONENT  ORIENTATION 


D.X  Introduction 

Given  the  body  (or  vehicle)  frame,  defined  by  x^,  Yfa,  Z^,  the 
(nominal)  inertial  component  frames  are  defined  such  that  the  input 
axes  of  one  inertial  component  case  of  each  kind  lies  along  (or  paral- 
lel to)  one  of  the  body  axes.  Hence  the  X-gyro  or  accelerometer  input 
axis  is  nominally  aligned  parallel  to  the  X-body  axis,  and  similarly 
for  the  other  inertial  components.  If  the  further  (and  usual)  con- 
straint - that  the  remaining  two  axes  of  each  inertial  component  be 
aligned  parallel  to  the  remaining  two  body  axes -is  applied,  then 
there  are  four  possible,  nominal  orientations  of  each  single-degree -of- 
frredow  inertial  component.  The  particular  orientations  chosen  are 
generally  predicated  on  an  attempt  to  minimize  the  degradation  in 
component  performance  in  the  anticipated  vehicle  environment,  bearing 

in  mind  that  calibration  is  usually  performed  in  a one  g field.  To 

2 

take  one  simple  example t if  a gyro  has  a large,  unstable,  g or  g de- 
pendent, torque  coefficient  which  is  forced  by  acceleration  (or  accel- 
eration squared)  spin  reference  axis  (SRA)  and  the  gyro  is 

to  be  used  in  an  aircraft  inertial  navigator,  then  the  first  choice  for 
the  orientation  of  the  gyro  SRA  would  be  along  the  lateral  or  Y-axis  of 
the  aircraft  since  there  is  normally  negligible  accleration  along  this 
axis.  The  second  choice  would  be  along  the  longitudinal  or  X-axis  of 
the  aircraft  where  there  is  minimal  or  only  transient  acceleration, 
the  final  choice  of  inertial  component  orientation  will  bo  a compromise 
based  on  several  similar  considerations. 
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The  default  orientation  of  the  inertial  components  appearing  in 
the  simulation  is  illustrated  in  Figure  D-l,  where  one  inertial  compo- 
nent of  each  kind  has  its  input  axis  parallel  to  one  positive  body  axis. 


D.2  X-inertial  Component  Orientation 

The  gyro  or  accelerometer  orientation  with  respect  to  the  body 
frame  is  specified  by  a rotation  matrix,  QGBX  or  QABX,  in  the  program 
mnemonics,  where  the  default  value  is: 

QGBX  or  QABX  - x(|) 

In  general,  the  nominal  orientations  for  the  X-gyro  or  accelerom- 
eter are: 

QGBX  or  QABX  « x(~],  * ° 0,  1,  2,  3 

If  one  wishes  to  specify  infinitesimal  misalignments  of  the  input  axis 
(XA)  of  the  particular  X inertial  component,  Xx>  with  respect  to  the 
body  frame,  these  may  be  represented  by  the  small  angles,  3xor  Yx» 
about  the  V and  2-body  axes.  The  resulting  transformation  from  the  body 
frame  to  the  slightly  misaligned  X-inertial  component  frame  is  given  by: 

QGBX  or  QABX  » x(^)ye(0x)  2e<Yx) 


10  0 
0 0 1 
0 10 


where  Y (0  ) and  Z (Yu)  denote  infinitesimal  rotations  about  the  Y and  Z 

EX  l X 

body  axes  respectively  and  are 


W 6V«V 


8 «i 

X 


1 o -8X 
0 10 
8,  o 1 
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<- 


and 


W - 2<V 


Y «1 
x 


0 

0 

1 


Since  $ and  Yv  are  small — on  the  order  of  a milliradian  or  leas— -their 
products  may  be  neglected  for  the  anticipated  applications,  and  the 
product  of  the  two  may  be  considered  to  be  multipiicatively  commuta- 
tive, i.e. , 


V.(6  ) 2^(Y  ) ^ 2^.<Y  ) <B  ) » 

e x e y e 'y  e x 


1 

“Y> 

8. 


Y -8 
x x 


1 

0 


0 
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The  general  transformation  from  the  body  frame  *-o  the  X inertial 
component  frame  is; 

“ 1 


QADX  or  gtibX  » 


-Y  C 
x 


(t)  Mt) 
♦v(t)  *At) 


c(t) 

«(t) 


I 


I 


-e» 

=(t) 

o(f) 


For  the  default  value  (k*i)  of  the  X inertial  component  orienta- 


tion 


QAUX  or  0<SSX 
(k  «*  1) 


1 

0 


0 

1 


of  course,  8 and  y may  be  different  for  the  X-gyro  than  for  the  X- 
x x 

accelerometer.  Misalignment  about  the  input  axis,  which  would  be  de- 
noted by  ax  in  this  case,  is  never  considered  in  practice. 
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D.3  Y-inertial  Component  Orientation 


The  V-gyro  or  accelerometer  orientation  with  respect  to  the  body 
frame  is  specified  by  a rotation  matrix,  QGBY  or  QABY,  in  the  program 
mnemonics,  where  the  default  value  is: 


QGBY  or  QABY  » X(ff) 


0 10 

il)  - 1 ° 0 

0 0-1 


In  general,  the  nominal  orientation  for  the  Y-gyro  or  accelerom- 
eter is: 

QGBY  or  QABY  ■ X^y)  k = 0,  1,  2,  3 

Infinitesimal  misalignments  of  the  XA  of  the  particular  Y-inertial 
component,  Xy,  with  respect  to  the  body  frame,  may  be  represented  by 
small  angles,  Yy  and  Cty,  about  che  2 and  X-body  axes  respectively.  The 
resulting  transformation  from  the  body  frame  to  the  slightly  misaligned 
Y-inertial  component  frame  is: 

QGBY  or  QABY  « x(~)z(f  )z£  (Yy)  X£Cay) 


where  z (y  ) and  Xt.  (o  ) denote  infinitesimal  rotations  about  the  Z and 

t y t y 

X-body  axes  respectively  and  are  s , 
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The  general  transformation  from  the  body  frame  to  the  Y -inertial 
component  frame  is: 


For  the  default  value  (k  » 2)  of  the  Y-inertial  component  orienta 

tion 


Finally*  the  2-gyro  or  accelerometer  orientation  with  respect  to 
the  body  frame  is  specified  by  a rotation  matrix,  QSB2  or  QAB2,  in  the 
program  mnemonics,  whore  the  default  value  isi 

QGB2  or  QAB2  « 2 (-  |)  y(-  |) 


in  general,  the  nominal  orientations  for  the  2-gyro  or  accelerom- 
eter are: 

QGBZ  OS  QAS2  » x[~r)  y(-  ~)  , k ■ 0,  1,  2,  3 


Infinitesimal  misa] ignments  of  the  IA  of  the  particular  Z-inertial 
component,  lz,  with  respect  to  the  body  frame  may  be  represented  by 
the  small  angles;  a and  $ , about  the  X and  Y body  axes  respectively. 

Z 2 

The  resulting  transformation  from  the  body  frame  to  the  slightly  mis- 
aligned Z-inertial  component  frame  is: 


QGBZ  or  fiUBZ  = x(f)v  -(f)*e«y 


where,  as  before 


W V6*>  68  V6z>  W 0 


1 

0 


0 

1 


-0, 


~BZ  -a. 


Hm  general  transformation  from  the  body  frame  to  the  Z-inortial 
component  frame  is: 
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D.5  Insertion  of  Inertial  Component  Orientation 


Each  element  of  each  inertial  component  alignment  matrix  must  be 
input  in  both  the  inertial  component  module  and  the  corresponding 
condensation  module.  It  should  be  reiterated  that  elements  of  these 
matrices  are  entered  row-wise,  not  columnwise;  for  instance,  if  QABZ, 
using  the  default  nominal  orientation  and  the  infinitesimal  misalign- 
ments, were  being  input  the  sequence  would  be: 


3 , , i,  l,  o,  -3  , 0,  l,  a . 

z z z z 


D.6  Component  Misalignments  and  Nonorthogonalities 

In  each  of  the  gyro  and  accelerometer  compensation  modules,  there 
is  a provision  for  an  additional  misalignment  matrix,  QMIS.  This 
matrix,  which  must  be  entered  by  the  user  as  a nine  element  vector,  is 
used  to  transform  the  compensated  AVas  and  A6gs  from  the  accelero- 
meter or  gyro  input  axes  to  the  body  frame.  (One  should  note  that 
QBAX  in  the  accelerometer  module  is  the  same  as  QABX  in  the  accelero- 
meter compensation  module  except  for  misalignment  errors,  i.e.  one  is 
not  supposed  to  be  the  inverse  of  the  other) . The  nine  elements  of 
QMIS  in  the  accelerometer  compensation  module  are  formed  from  the 
elements  of  the  first  rows  of  each  of  QABX,  QABY,  and  QABZ  with  the 
signs  of  the  infinitesimal  misalignments  reversed,  i.e., 


1 -y  3 
x x 


QMIS  A 


-6,  a, 

. z z 


m 

qabx^ 

-QABX2 

-QABX3 

-QABY^ 

Qaby„ 

d. 

-QABY, 

_ -QAB21 

-QABZ2 

QABZ3  _ 

Of  course,  the  diagonal  elements  of  yMIS  are  all  unity.  QMIS  is  the 
"first  order  inverse"  of  the  nonorthogonal  transformation  from  the  body 
frame  to  the  input  axes  of  the  accelerometer  triad.  Shortening  QMIS  to 
Q,  the  following  identity  holds: 

Q \ (Q  + Q)  + ’2  (0  - Q) 


-Y  + Y 
x y 


, -H  t B 

L X 


-Y  + v 
x y 


-a  + a 
y z 


-B  + 6 
z x 


*n  + n 
Y : 


Symmetric 


(Yx  4 V 


Hz  + 
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Y t Y 
x y 


,-(B  + B ) 

L Z X 


0 -(a  + a ) 

Y z 


« + o 0 

Y z 


Skew 

Symmetric 


The  off  diagonal  element  of  the  symmetric  portion  of  Q represent  the 
component  non  orthogonalities.  The  skew  symmetric  portion  of  Q repre- 
sents the  rotational  misalignment  of  the  orthoqonalized  triads. 

It  should  be  noted  that  the  a's,  P's,  and  y's.  do  not  necessarily 
imply  the  existence  of  component  alignment;  only  when  the  values  of  Yy 
and  p^  in  yABX,  for  instance,  do  not  coincide  with  the  corresponding 
values  in  QBAX  do  actual  (first  order)  errors  in  alignment  of  the  X 
accoloromoter  propagate  into  systom  orrors. 


IV) 


The  form  of  the  component  misalignments  used  throughout  the  INSS 
do  imply  the  existence  of  independent  means  of  establishing  the  'body* 
frame.  Tnis  in  turn  implies  the  existence  of  optical  cubes,  one  of 
which  defines  the  "body"  frame  fixed  in  the  inertial  component  cluster. 

As  it  stands,  it  is  evident  that  the  user  must  exercise  great 
care  in  setting  up  a dosired  set  of  nyro  and  accelerometer  nonortho- 
gonalities  and  rotational  misalignments. 
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